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This  grant  consisted  of  work  divided  into  two  parts.  The  first  part,  “Ocean 
Modeling,”  was  performed  off-campus.  The  lead  researcher  was  Dr.  Germana 
Peggion  of  the  Department  of  Physics  at  UNO.  Her  colleagues  for  this  research 
were  employees  of  the  Oceanography  Division  of  the  Naval  Research  Laboratory 
at  Stennis  Space  Center  (NRL-SSC).  The  main  objective  of  the  ocean  modeling 
portion  of  this  grant  was  to  collaborate  with  and  assist  NRL-SSC  in  implementing, 
evaluating,  and  applying  ocean  forecasting  systems  in  support  of  naval 
applications. 

The  second  part  of  the  grant  (Coastal  Modeling)  consisted  of  calculations 
concerning  storm  surge  simulations  over  south  Louisiana  and  was  performed  on 
campus.  The  lead  investigator  was  Dr.  Martin  Guillot  of  the  Department  of 
Mechanical  Engineering  at  UNO.  He  was  assisted  in  this  work  by  Dr.  Ioannis 
Georgiou  of  the  Pontchartrain  Institute  for  Environmental  Sciences  at  UNO  and 
Dr.  Alex  McCorquodale  of  the  Department  of  Civil  and  Environmental 
Engineering  at  UNO.  The  objective  was  to  compute  water  surface  elevations  over 
southern  Louisiana  due  to  storm  surge  produced  by  the  updated  “Standard  Project 
Hurricane”  (SPH).  The  SPH  defines  a  methodology  for  producing  hurricane  path, 
forward  speed  and  wind  field  scenarios  that  are  used  as  design  criteria  to  guide  the 
design  heights  of  the  levee  system  in  southeast  Louisiana. 

Dr.  George  Ioup  and  Dr.  Juliette  Ioup  of  the  Department  of  Physics  at  UNO 
administered  the  grant. 

There  are  five  parts  to  this  report.  The  first  is  this  overview.  The  second  is  a 
summary  of  achievements  for  the  Ocean  Modeling  by  Dr.  Peggion  organized  by 
Task  and  titled  Ocean  Modeling.  The  third  is  a  summary  of  the  achievements  for 
the  Coastal  Modeling  by  Dr.  Guillot  and  Dr.  Georgiou  titled  Standard  Project 

20090401083 


Hurricane  Update:  ADCIRC  Storm  Simulations  Over  Southeast  Louisiana.  The 
fourth  and  fifth  are  two  manuscripts  to  be  published  in  the  Journal  of  Marine 
Systems,  which  summarize  some  of  the  work  by  Dr.  Peggion  and  her  colleagues, 
Super-ensemble  Forecasts  and  Resulting  Acoustic  Sensitivities  in  Shallow  Waters, 
and  A  Note  on  Ncom  Temperature  Forecast  Error  Calibration  Using  the  Ensemble 
Transform.  The  latter  papers  were  also  supported  by  a  funding  from  NRL-SSC 
through  NASA-Stennis  and  the  Louisiana  Board  of  Regents  to  UNO  as  well  as  by 
a  follow-on  NRL  cooperative  agreement  with  UNO,  N00173-07-2-C901. 


STANDARD  PROJECT  HURRICANE  UPDATE:  ADCIRC 
STORM  SIMULATIONS  OVER  SOUTHEAST  LOUISIANA 

By 

Dr.  Martin  J.  Guillot 
Dr.  Ioannis  Georgiou 
University  of  New  Orleans 
New  Orleans,  LA 


Objective 

The  purpose  of  the  current  effort  is  to  compute  storm  surges  produced  by  the 
standard  project  hurricane  (SPH)  using  the  ADCIRC  storm  surge  model  and  to  compare 
surges  resulting  from  indices  defined  in  a  1959  National  Weather  Bureau  report  with 
indices  defined  as  part  of  the  SPH  reanalysis  after  the  2005  hurricane  season. 

Standard  Project  Hurricane 

The  SPH  is  one  of  the  design  criteria  the  U.S.  Army  Corps  of  Engineers  uses  for 
hurricane  protection  projects  along  the  east  and  gulf  coasts  of  the  United  States.  The  SPH 
is  a  hypothetical  hurricane  based  on  historical  data  and  observations  of  hurricanes  that 
have  occurred  the  Atlantic  basin  along  the  east  and  gulf  coasts  and  was  originally  defined 
in  1959  as  part  of  the  National  Hurricane  Research  Project  in  the  U.S.  Weather  Bureau 
Report  No.  33  [1]  (subsequently  referred  to  as  NHRP  33)  using  data  from  storms  that 
occurred  during  the  period  1900-1956.  That  report  defined  the  SPH  as  “...the  most  severe 
storm  that  is  considered  reasonably  characteristic  of  a  region  in  which  the  basin  is 
located”.  The  U.S.  Weather  Bureau  and  the  U.S.  Army  Corps  of  Engineers  (USACE) 
jointly  derived  the  specifications,  criteria  and  procedures  for  computing  the  SPH  defined 
in  NHRP  33.  The  U.S  east  and  gulf  coasts  were  divided  into  zones  of  approximately 
equal  area  and  SPH  indices  were  defined  for  each  zone.  The  three  zones  defined  on  the 
gulf  coast  are  shown  in  Figure  1  (taken  from  NHRP  33).  The  primary  indices  used  to 
define  the  SPH  within  each  zone  are:  central  pressure  index  (CPI),  maximum  30  ft  (10  m) 
over  water  winds,  radius  of  maximum  winds  and  forward  speed.  The  New  Orleans  area  is 
located  in  zone  B. 

The  SPH  has  undergone  reanalysis  several  times  since  originally  defined  in 
NHRP  33.  Reanalysis  has  included  redefining  both  the  SPH  indices  based  on  more  recent 
(after  1956)  data  and  refining  the  methodologies  used  to  compute  the  wind  and  pressure 
fields.  In  1979  a  reanalysis  of  the  SPH  indices  based  on  storms  through  1975  was 
published  in  NOAA  Technical  Report  23  [2]  (subsequently  referred  to  as  NWS  23).  That 
document  redefined  the  SPH  as  “a  steady  state  hurricane  having  a  severe  combination  of 
values  of  meteorological  parameters  that  will  give  high  sustained  wind  speeds  reasonably 
characteristic  of  a  given  region”,  and  also  revised  the  methodologies  used  for  computing 
the  wind  and  pressure  fields  from  those  used  in  1959  NRRP  33  report.  In  1996, 

Thompson  and  Cardone  [5]  developed  a  model  for  generating  tropical  cyclones  based  on 
the  planetary  boundary  layer  approach.  This  approach  was  incorporated  into  the  Ocean 


Weather,  Inc  (OWI)  meso-scale  vortex  numerical  model  for  specification  of  surface  wind 
and  pressure  fields  inside  tropical  cyclones  based  on  specifying  the  appropriate  indices. 
After  the  2005  hurricane  season,  the  SPH  indices  were  reanalyzed  again  to  include  all 
data  from  1851  through  the  2005  hurricane  season.  Levinson  [3]  presents  the  results  of 
the  SPH  reanalysis  and  compares  the  1979  indices  defined  in  NWS  23  with  the  new 
indices  defined  as  part  of  the  SPH  reanalysis  after  the  2005  hurricane  season. 

The  current  study  focuses  on  comparing  storm  surges  predicted  by  the  ADCIRC 
storm  surge  model  using  the  SPH  indices  defined  in  NHRS  33  (1959)  and  the  SPH 
indices  defined  in  the  reanalysis  after  the  2005  hurricane  season.  For  the  remainder  of 
this  report,  the  SPH  indices  as  defined  in  NHRS  33  will  be  referred  to  as  the  “old”  SPH 
and  the  SPH  indices  defined  as  part  of  the  SPH  reanalysis  after  the  2005  hurricane  season 
will  be  referred  to  as  the  “new”  SPH.  The  methodology  for  computing  the  wind  and 
pressure  fields  for  the  old  SPH  is  based  on  the  methods  in  the  1979  NWS  23  report  and 
the  methodology  for  computing  the  wind  and  pressure  fields  for  the  new  SPH  is  based  on 
OWI  tropical  cyclone  model  called  TC96.  Storm  surges  are  computed  using  the  shallow 
water  modeling  system  ADCIRC 

The  approximate  SPH  wind  field  parameters  used  in  this  study  are  presented  in 
Table  1. 


Table  I:  SPH  indices  used  to  define  the  new  and  old  SPH. 


Central 
Pressure  Index 
(mb) 

Radius  of 
Maximum  Winds 
(nautical  mile) 

Max  Wind 
Speed 
(mjph) 

Method  for 
computing  wind 
and  pressure 

New  SPH 

904.1 

11 

132.0 

OWI  TC96  model 

Old  SPH 

934.6 

30 

104.0 

NWS  23 

Figure  1:  Gulf  coast  zones  defined  in  Report  33. 


ADCIRC  Storm  Surge  Modeling  System 

The  ADCIRC  storm  surge  model  was  originally  developed  by  Leuttich  and 
Westerink  [4]  and  since  then  has  undergone  extensive  development  by  several 
researchers  and  organizations.  The  model  is  based  on  the  finite  element  method  and 
incorporates  the  generalized  wave  continuity  equation  (GWCE)  for  numerical  stability. 
ADCIRC  computes  water  surface  elevations  and  velocities  at  nodal  points  for  the  2-D 
depth  averaged  shallow  water  equations.  A  parallel  version  of  ADCIRC  has  been 
developed  to  run  on  several  parallel  architectures,  including  Linux  clusters.  For  this  study 
the  ADCIRC  model  was  run  on  the  64  node  Linux  cluster  at  the  University  of  New 
Orleans.  Version  46.52  is  used  in  conjunction  with  the  southeast  Louisiana  mesh 
sll5v3  2007  r09a.  The  code,  mesh,  control  files  as  well  as  wind  fields  were  provided  by 
the  U.S.  Army  Corps  of  Engineers,  New  Orleans  District.  The  wind  fields  used  to  force 
the  ADCIRC  model  and  are  discussed  in  more  detail  below. 

Modifications  to  Original  Mesh 

The  ADCIRC  mesh  consists  of  nodes,  elements  and  boundaries.  ADCIRC  has  the 
capability  to  model  several  types  of  boundary  conditions  including,  but  not  limited  to, 
elevation  (tidal),  inflow  (river),  outflow  and  weir.  Levees  are  modeled  using  weir 
boundaries.  Weirs  are  specified  at  the  given  levee  height  for  each  levee  using  weir  node 
pairs  and  when  the  water  surface  elevation  exceeds  the  elevation  specified  at  the  weir 
node  pair,  weir  equations  are  used  to  compute  flow  over  that  levee  to  simulate  levee 
overtopping.  For  the  SPH  study,  levee  overtopping  was  prevented  in  the  Lake 
Pontchartrain  and  vicinity,  including  the  west  bank  by  specifying  the  elevation  of  the 
node  weir  pairs  of  the  flood  protection  levees  to  be  25  m.  The  purpose  is  to  help 
determine  levee  heights  that  would  be  required  to  protect  the  New  Orleans  area  and  west 
bank  from  an  SPH  event.  The  specific  levees  raised  are  shown  in  Figure  2. 


Figure  2;  Node  weir  pairs  that  simulate  flood  protection  levees  in  and  around  New  Orleans  that  were 
raised  to  25  m. 


SPH  Tracks  and  Wind  Fields 

The  SPH  study  consists  of  three  primary  tracks  labeled  A,  C  and  F.  The  tracks  are 
based  on  COE  experience  and  historical  storm  data,  and  are  shown  in  Figure  3  along  with 
the  forward  translational  speed  along  each  track.  The  old  and  new  SPH  wind  field 
contours  are  shown  Figure  4  and  Figure  5,  respectively.  Several  things  are  noted  about 
the  wind  fields.  While  the  new  SPH  is  a  stronger  storm  in  terms  of  maximum  wind  speed, 
it  is  also  a  much  smaller  storm,  with  radius  of  maximum  winds  almost  one  third  of  the 
old  SPH.  Additionally,  as  can  be  seen  from  comparing  the  two  figures,  the  wind  field 
extends  farther  for  the  old  SPH  than  the  new  SPH.  For  example,  for  the  storm  positions 
shown,  the  wind  speed  in  the  Mississippi  sound  near  the  entrance  to  Lake  Borgne 


Standard  Project  H  urn  cane  Tracks 

Track  A  -  OLD,  Forward  Translation  6  knota 
— —  Track  A  -  NEW,  F orward  Translation  6  knots 

Track  C-  OLD,  Forward  Tranalation  6  krtota 
■  ■  Track  C  -  NEW,  Forward  Translation  6  knota 

-  Track  F  -  OLD,  Forward  Tranalation  1 1  knota 

. .  Track  F -NEW,  Forward  Translation  1 1  knots 


Figure  3:  SPH  Storm  Tracks  and  forward  translational  speed. 


Figure  4:  Old  SPH  wind  field 


Figure  5:  New  SPH  wind  field 


is  on  the  order  of  70-80  mph  for  the  old  SPH,  whereas  for  the  new  SPH,  the  wind  speed 
has  dropped  to  40-50  mph.  The  larger  size  of  the  old  SPH  can  have  a  significant  effect  on 
the  computed  storm  surge,  as  will  be  shown  in  the  results  section. 

Results 

The  results  of  the  study  are  presented  in  a  series  of  contour  plots  of  maximum 
water  surface  elevation,  elevation  difference  and  water  surface  elevation  histories  at 
selected  locations  within  the  Lake  Pontchartrain  and  vicinity  study  area. 

All  simulations  were  begun  30  hours  before  land  fall  with  wind  ramping  to  full 
strength  during  the  first  6  hours,  so  that  full  force  winds  began  24  hours  prior  land  fall. 
Prior  to  performing  the  SPH  simulations,  ADCIRC  was  run  without  winds  for  a  2  day 
simulation  with  inflow  boundary  conditions  specified  on  the  Atchafalaya  and  Mississippi 
rivers  to  provide  correct  initial  water  surface  elevations  for  the  SPH  simulations  These 
results  were  saved  and  the  SPH  simulations  were  hot  started  from  these  initial  conditions. 
No  tidal  dynamics  were  included  in  this  study. 

The  maximum  water  surface  elevation  contours  for  tracks  A,  C  and  F  are  shown  in  Figure  6 
through  Figure  8  for  the  (a)  old  and  (b)  new  SPH  respectively.  Figure  6  shows  that  for  track  A,  the 
surge  produced  along  the  south  shore  of  Lake  Pontchartrain  is  higher  for  the  old  SPH  than  for  the 
new  SPH.  This  is  due  to  the  larger  area  of  the  wind  field  for  the  old  SPH.  Near  the  MRGO  and 
Plaquemines  Parish  on  the  east  side  of  the  Mississippi  river,  the  storm  surges  for  the  old  and  new 
SPH  are  comparable.  However,  the  higher  surge  extends  farther  into  the  Mississippi  sound  for  the 
new  SPH.  This  is  expected,  since  the  eye  passes  over  this  area  and  the  maximum  winds  near  the  eye 
wall  are  higher  for  the  new  SPH.  For  track  C  it  is  clearly  seen  by  comparing 

Figure  7  (a)  and  (b)  that  the  old  SPH  produces  a  higher  storm  surge  for  most  of  the  area 
of  interest,  except  on  the  south  shore  of  lake  Pontchartrain  around  New  Orleans  east. 
Figure  8  (a)  and  (b)  show  that  the  storm  surge  produced  by  track  F  is  substantially  lower 
than  for  track  A  and  C.  Figure  9  (a)  and  (b)  shows  the  maximum  water  surface  elevations 


produced  by  considering  all  tracks  for  the  old  and  new  SPH,  respectively.  Figure  10 
through  Figure  1 3  show  the  differences  in  computed  water  surface  elevation  (new  -  old) 
for  tracks  A,  C,  F  and  considering  the  maximum  of  all  tracks,  respectively.  The 
difference  contour  plots  more  clearly  demonstrate  the  differences  in  computed  water 
surface  elevation  between  the  old  and  new  SPH. 

Next,  water  surface  elevation  histories  were  computed  at  selected  points  in  the 
flood  protection  system  in  the  New  Orleans  and  surrounding  areas.  The  elevation 
recording  stations  are  shown  in  Figure  14  and  are  grouped  for  plotting  purposes  into  four 
sections  along  the  flood  protection  system.  The  maximum  water  surface  elevations  for  at 
each  of  the  22  stations  are  shown  in  Table  ()  for  each  storm  and  the  maximum  of  all 
storms. 

Figure  15  shows  that  along  the  south  shore  of  Lake  Pontchartrain,  the  old  SPH 
produces  a  maximum  storm  surge  of  approximately  5.5  meters  whereas  the  new  SPH 
produces  a  maximum  storm  surge  between  4  and  5  meters.  For  both  storms  the  storm 
surge  is  relatively  constant  at  stations  2,  3  an  4,  but  is  higher  at  stations  1  and  5.  Also,  the 
storm  surge  peaks  approximately  2  hours  later  for  the  old  SPH.  Figure  16  indicates  that 
the  maximum  storm  surge  in  the  New  Orleans  east  section  occurs  at  station  10  and  is 
approximately  7  meters.  In  the  MRGO/IHNC  section,  Figure  17  indicates  that  a 
maximum  storm  surge  of  8  meters  occurs  at  stations  13-16  and  is  approximately  equal  for 
both  storms.  No  appreciable  flooding  occurred  at  the  west  bank  stations  for  track  A  and 
so  those  stations  are  not  plotted. 

Figure  18  shows  that  track  C  produces  significantly  lower  peak  surge  along  the 
south  shore  of  Lake  Pontchartrain  than  track  A,  with  a  maximum  surge  occurring  at 
station  1  of  approximately  4.5  meters  for  the  old  SPH  and  just  under  4  meters  for  the  new 
SPH.  In  New  Orleans  east.  Figure  19  shows  that  the  maximum  storm  surge  occurs  at 
station  10  and  is  6.4  meters  for  the  old  SPH  and  4.4  meters  for  the  new  SPH.  Figure  20 
shows  maximum  surge  occurring  at  station  1 1  of  approximately  8  meters  for  the  old  SPH 
and  6.4  meters  for  the  new  SPH.  Stations  13-16  show  a  maximum  surge  of  7.8  meters  for 
the  old  SPH  and  5.8  for  the  new  SPH.  On  the  west  bank,  Figure  21  indicates  that  the  old 
SPH  produces  significantly  higher  storm  surge  at  stations  18  and  19  with  at  maximum 
surge  of  7.8  meters  for  the  old  SPH  and  a  maximum  surge  of  5.0  meters  for  the  new 
SPH.  Unfortunately,  stations  20-22  indicate  that  the  surge  has  not  yet  peaked  for  the  old 
SPH,  indicating  that  the  simulation  ended  before  the  maximum  surge  was  captured.  The 
new  SPH  produces  a  maximum  surge  of  less  than  4  meters  at  those  locations. 

Elevation  histories  produced  by  Track  F  are  shown  in  Figure  22  through  Figure 
25.  In  contrast  to  Tracks  A  and  C,  the  new  SPH  produces  higher  storm  surges  for  Track  F 
rather  than  the  old  SPH  at  the  selected  locations.  This  is  not  surprising,  since  the  eye  of 
the  storm  passes  over  the  city  and  the  new  SPH,  although  smaller  has  stronger  winds  near 
the  eye.  It  is  seen  from  Figure  22  that  Track  F  produces  significant  storm  surge  on  the 
south  shore  of  Lake  Pontchartrain  and,  in  fact,  produces  higher  storm  surges  there  than 
tracks  A  and  C.  At  the  other  locations  the  storm  surge  is  still  significant,  but  not  as  high 
as  Track  A  or  C. 

The  maximum  water  surface  elevations  at  each  recording  station  for  each  storm 
track  are  presented  in  Table  2  as  well  as  the  overall  maximum  for  both  the  old  and  new 
SPH. 


Conclusions 

The  study  produced  some  unexpected  results.  Initially  it  was  thought  that  the 
stronger  storm  defined  by  the  new  SPH  would  produce  higher  storm  surges,  but  that  was 
not  always  the  case.  It  is  clear  from  the  maximum  water  surface  elevation  contours  and 
water  surface  elevation  histories  that  the  larger  storm  defined  by  the  old  SPH  indices  has 
significant  influence  on  the  resulting  peak  storm  surges  at  most  locations  of  interest. 


Figure  6:  Computed  maximum  water  surface  elevation  SPH  Track  A,  (a)  old  SPH,  (b)  new  SPII 


Figure  7:  Computed  maximum  water  surface  elevation  SPH  Track  C,  (a)  old  SPH,  (b)  new  SPH 


Figure  9:  Computed  maximum  water  surface  elevation,  locus  of  Tracks  A,  C  and  F,  (a)  old  SPH,  (b) 
new  SPH 


Figure  10:  Difference  in  computed  maximum  water  surface  elevation  (new-old),  Track  A 


Figure  11:  Difference  in  computed  maximum  water  surface  elevation  (new-old),  Track  C. 


Figure  12:  Difference  in  computed  maximum  water  surface  elevation  (new-old),  Track  F. 


Figure  13:  Difference  in  computed  maximum  water  surface  elevation  (new-old),  locus  of  tracks  A,C 
and  F. 
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Figure  14:  Elevation  recording  station  locations  grouped  into  sections. 
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Figure  15:  Track  A  SPH  Elevation  histories  at  selected  locations  along  south  shore  of  Lake 
Pontchartrain  in  St.  Charles  and  Jefferson  Parish  and  New  Orleans  Lakefront. 
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Figure  16:  Track  A  SPH  Elevation  histories  at  selected  locations  along  south  shore  of  Lake 
Pontchartrain  in  New  Orleans  East. 
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Track  A  SPH  Elevation  histories  at  selected  locations  along  MRGO,  industrial  canal,  and 
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Figure  18:  Track  C  SPH  Elevation  histories  at  selected  locations  along  south  shore  of  Lake 
Pontchartrain  in  St.  Charles  and  Jefferson  Parish  and  New  Orleans  Lakefront. 
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Figure  19:  Track  C  SPH  Elevation  histories  at  selected  locations  along  south  shore  of  Lake 
Pontchartrain  in  New  Orleans  East. 
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Figure  20:  Track  C  SPH  Elevation  histories  at  selected  locations  along  MRGO,  industrial  canal,  and 
GIWW. 
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Figure  21:  Track  C  SPH  Elevation  histories  at  selected  locations  on  West  Bank. 
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Figure  22:  Track  F  SPH  Elevation  histories  at  selected  locations  along  south  shore  of  Lake 
Pontchartrain  in  St.  Charles  and  Jefferson  Parish  and  New  Orleans  Lakefront. 
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Figure  23:  Track  F  SPH  Elevation  histories  at  selected  locations  along  south  shore  of  Lake 
Pontchartrain  in  New  Orleans  East. 
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Figure  24:  Track  F  SPH  Elevation  histories  at  selected  locations  along  MRGO,  industrial  canal,  and 
GIWW. 
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Figure  25:  Track  F  SPH  Elevation  histories  at  selected  locations  on  West  Bank. 


Table  2:  Maximum  recorded  water  surface  elevations  at  selected  recording  stations. 


Maximum  Water  Surface  Elevation  at  Recording  Stations 

Track  A 

Track  C 

Track  F 

Maximum 

Station 

Old 

New 

Old 

New 

Old 

New 

Old 

New 

1 

5.4 

4.3 

4.7 

4.1 

3.7 

5.7 

5.4 

5.7 

2 

5.5 

4.2 

4.1 

3.7 

3.7 

5.7 

5.5 

5.7 

3 

4.8 

3.5 

3.7 

2.9 

3.0 

4.9 

4.8 

4.9 

4 

4.8 

3.4 

3.0 

2.4 

2.7 

4.7 

4.8 

4.7 

5 

5.2 

3.5 

2.1 

1.7 

1.7 

4.1 

5.2 

4.1 

6 

5.3 

3.5 

1.0 

0.8 

1.6 

2.8 

5.3 

3.5 

7 

5.0 

3.3 

2.0 

2.0 

1.5 

1.5 

5.0 

3.3 

8 

4.7 

3.5 

2.6 

2.0 

1.8 

1.7 

4.7 

3.5 

9 

6.8 

6.7 

5.3 

3.5 

2.9 

3.7 

6.8 

6.8 

10 

9.1 

9.1 

6.8 

5.0 

4.6 

5.2 

9.1 

9.1 

11 

8.0 

8.0 

8.0 

7.0 

2.7 

3.2 

8.0 

8.0 

12 

6.5 

6.7 

6.1 

4.8 

3.7 

5.1 

6.5 

6.7 

13 

8.0 

8.5 

7.5 

5.5 

4.4 

6.5 

8.0 

8.5 

14 

7.9 

8.2 

7.4 

5.5 

4.1 

6.7 

7.9 

8.2 

15 

7.8 

8.1 

7.2 

5.4 

4.0 

6.5 

7.8 

8.1 

16 

7.7 

8.0 

7.1 

5.4 

3.9 

6.5 

7.7 

8.0 

17 

6.0 

6.2 

5.5 

4.3 

2.7 

5.0 

6.0 

6.2 

18 

0.0 

0.0 

7.5 

5.0 

1.8 

1.8 

7.5 

5.0 

19 

0.0 

0.0 

7.3 

4.8 

1.2 

1.2 

7.3 

4.8 

20 

0.0 

0.0 

4.5 

3.3 

0.7 

1.1 

4.5 

3.3 

21 

0.0 

0.0 

4.5 

3.6 

1.7 

1.5 

4.5 

3.6 

22 

0.0 

0.0 

3.0 

3.6 

0.0 

0.0 

3.0 

3.6 
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Ocean  Modeling 


Task  1.  Model  Hindcasting  in  Smart  Climatology 

Collaboration  with:  James  Dykes1,  Lucy  Fitzgerald  Smedstad1 

Research  Objectives:  This  project  aims  to  develop  concepts  and  techniques  that  clearly  define 
smart  climatology  for  strategic  planning  for  ASW,  SPECOPS  and  MIW,  enabling  superior 
analysis  of  environmental  variability  to  support  tactical  decision  planning.  Ultimately,  smart 
climatology  is  to  take  into  account  the  effects  on  strategic  planning  of  tactical  extremes  in  ocean 
and  atmospheric  conditions  caused  by  the  large-scale  climatic  variations.  A  demonstration  of 
this  project  is  planned  to  provide  the  guidance  for  future  potential  transition  of  an  end-to-end 
capability  to  the  war  fighter.  The  portion  of  this  project  requiring  HPC  resources  involves 
running  (MetOc)  models  for  long  historical  periods  with  the  goal  of  providing  quick-turn- 
around  results  on  demand. 

Methodology:  MSRC  resources  were  utilized  in  generating  the  data  and  information  based  on 
running  atmospheric  and  oceanographic  models  over  a  long  period  of  time  in  the  past,  also 
known  as  hindcasts.  A  high  resolution  tactical  scale  climatology  dataset  required  for  knowledge 
extraction  was  generated  by  an  air/ocean/wave  coupled  system,  which  has  been  constrained  by 
relevant  large-scale  climatic  variations.  The  system  components  include  COAMPS  , 
WAVEWATCH  III,  NCOM,  and  SWAN  for  creating  strategic  and  tactical  climatologies  in 
data-sparse  and  data-void  areas,  creating  a  three-dimensional  depiction  of  the  atmosphere  and 
the  ocean  over  a  three-year  period  (1997  through  1999).  This  period  is  limited  in  time  to  cover 
the  anomalous  events  of  an  extreme  El  Nino  and  La  Nina  for  demonstration,  and  is  expected  to 
expand  in  later  work.  Certain  parameters  are  extracted  depending  on  the  mission  scenario. 
Ultimately,  all  the  models  will  be  closely  coupled  under  ESMF,  but  for  now  they  were  run 
separately.  The  ocean  models  used  forcing  provided  by  either  NCEP/NCAR  Reanalysis  or 
COAMPS*  run  at  NRL  Monterey.  Global  NCOM  output  provided  the  boundary  conditions  for 
the  regional  NCOM. 

UNO  contribution:  NCOM  in  a  2  coupled  nest  configuration  were  run  from  midway  in  1997 
through  1999.  All  the  resulting  model  output  files  including  the  complete  atmospheric  model 
outputs  were  stored  on  the  Sun-Fire- 15000  (vincent)  at  NAVO  MRSC  server  occupying  about 
15  terabytes  total.  This  server  provided  a  convenient  means  for  data  sharing  amongst  team 
members.  The  processed  output  was  passed  on  to  NRL-7440  to  be  used  in  pattern  analysis 
procedures  resulting  in  information  that  will  allow  to  examine  and  validate  the  types  of  data 
and  statistics  that  may  impact  strategic  planning. 


Task  2.  Relocatable  NCOM 

Collaboration  with:  C.  Rowley,  R.  Allard,  E.  Coelho 

Research  Objectives:  The  main  purpose  of  this  project  is  to  develop  and  evaluate  a  real  time 
ocean  prediction  system  developed  at  the  Naval  Research  Laboratory  (NRL)  in  support  of  naval 
operations.  The  system  is  portable  on  several  computer  platforms  and  operating  systems,  and 


rapidly  relocatable.  Analysis  and  prediction  are  available  for  any  part  of  the  world  usually 
within  a  few  hours  from  the  request,  making  it  a  particularly  useful  system  in  emergency 
situations.  The  major  challenge  is  to  offer  a  default  set  of  parameters  that  can  provide  accurate 
solutions  for  any  given  configuration,  yet  allowing  the  flexibility  of  tuning  and  calibrating  for  a 
given  domain  configuration.  Since  is  it  unrealistic  to  assume  data  are  available  at  the  spatial 
and  temporal  resolutions  necessary  for  specification  of  the  boundary  conditions,  the  system  has 
the  capability  of  multiple  1-way  nesting  from  basin-scale  to  regional  to  high-resolution  coastal 
domains. 

Methodology:  The  relocatable  system  was  evaluated  in  several  realtime  configurations  in 

support  of  NAVO  operations  and  other  research  joint  program.  While  the  exercises  at  NAVO 
have  the  support  of  allocated  computer  resources,  other  realtime  applications  are  sensibly 
constrained  by  the  computational  requirement  and  a  timely  deliver  of  the  solution 

UNO  contributions:  in  FY  2006-2007,  the  relocatable  system  was  evaluated  in  2  major 
realtime  exercises: 

1.  SW_06  (Shallow  Water  2006)  a  joint  experiment  with  NRL  and  other  academic 
institutions  (Rutgers  University  being  the  leader  organization)  off  the  New  Jersey  coast. 
Relo  NCOM  was  configured  in  3  nest  domains  with  a  horizontal  resolution  ranging 
from  2.5  to  .6  km.  the  inner  nest  was  designed  to  capture  the  internal  wave  activity  at 
the  shelf  break  and  provide  accurate  forecast  to  the  acoustic  group. 

2.  MREA_07  (Marine  Rapid  Environmental  Assessment  2007).  One  of  the  NRL 
contributions  to  the  exercise  was  to  provide,  in  realtime,  ocean  forecasts  in  support  of 
the  operations  at  sea.  The  NRL  prediction  system,  was  configured  with  3  nesting 
domains  at  resolutions  of  4,  2,  and  0.6  km.  Two  separate  inner  nests  were  configured 
for  the  BP  07  (Elba)  and  LASIE  (LaSpezia)  areas  of  operations,  respectively.  For  this 
application,  no  data  were  assimilated  in  realtime.  However,  a  small  (10  members) 
ensemble  of  free-runs  was  used  for  water  column  temperature  forecast  Root  Mean 
Square  (RMS)  error  prediction.  Ocean  forecast  arc  usually  the  final  component  of  a 
long  string  of  products  developed  at  several  different  centers:  a  delay  in  acquiring  one  of 
the  input  data,  the  classic  computer  breakdowns  (just  to  mention  a  few  issues)  may 
create  a  domino  effect  and  ultimately  a  late  delivery  of  the  forecast.  Preliminary 
model/data  comparison  and  new  simulations  in  a  pseudo  forecast  mode,  but  with 
different  model  parameters  (such  as  increased  vertical  resolution)  highlight  the  skills  and 
limits  of  the  default  configuration. 

Presentations: 

G.  Peggion:  How  to  develop  a  relocatable  prediction  system.  OGS  Trieste,  May  19  2006 
G.  Peggion:  A  Rapidly  Relocatable  Ocean  Prediction  System:  Congress  SIMAI  2006.  Baia  S. 
Samuele,  Italy, May,  25,2006 
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Lerici,  Italy.  Sept  25-27,  2007 
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A  NOTE  ON  NCOM  TEMPERATURE 


FORECAST  ERROR  CALIBRATION 
USING  THE  ENSEMBLE  TRANSFORM 


Emanuel  F.  Coelhol,2\  Germana  Peggion1  \  Clark  Rowley1,  Gregg  Jacobs1,  Richard 

Allard1  and  Elaina  Rodriguez3 

Naval  Research  Laboratory,  Stennis  Space  Center,  MS-USA 
University  of  Southern  Mississippi,  MS-USA 
^University  of  New  Orleans,  Louisiana,  LA-USA 


Abstract: 

During  the  MREA07  trial,  off  the  NW  coast  of  Italy  in  the  late  spring  and  summer  of  2007,  Navy 
Coastal  Ocean  Modeling  (NCOM)  multiple  nests  free  run  ensembles  were  made  available  in  real¬ 
time  for  the  LASIE07  and  BP07  events  and  a  fairly  complete  set  of  observations  were  collected 
inside  the  inner  nests  domains.  This  note  addresses  the  problem  of  predicting  NCOM  local 
unbiased  0-24  hours  forecast  errors  by  perturbing  a  limited  number  of  possible  error  sources 
through  Monte-Carlo  simulations,  without  local  data  assimilation.  It  discusses  preliminary  results 
using  the  Ensemble  Transform  (Bishop  and  Toth,  1999)  to  calibrate  the  ensemble  spread  by 
adjusting  its  characteristics  (spread-skill  relationship  and  magnitude)  to  an  observed  or  pre¬ 
estimated  error  field.  A  small  (10  members)  ensemble  of  free-runs  was  used  for  water  column 
temperature  forecast  Root  Mean  Square  (RMS)  error  prediction.  After  being  post-processed  they 
were  compared  with  observed  errors  and  those  estimated  using  time  variability  as  an  error  proxy. 
The  ensemble  runs  were  generated  through  atmospheric  forcing  perturbations  using  the  space- 
time  deformation  method  as  proposed  by  Xiandong,  et  al.,  (2007),  keeping  independent  initial 
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conditions.  Because  at  the  starting  time  all  runs  shared  the  same  IC,  the  ensemble  was  run  for 
roughly  two  weeks  for  spinning  up  and  then  used  during  the  following  10  days  for  data 
comparisons,  dunng  which  the  ensemble  spread  did  not  diverge  and  was  consistent  with  the 
observed  dynamics.  Comparisons  of  ensemble  spread  of  temperature  profiles  with  local  observed 
errors  and  time  variability  (assumed  as  an  error  proxy)  showed  they  were  consistent  through  this 
10  days  analysis  period,  with  performances  above  the  non-calibrated  ensemble  estimates  and 
time-variability  used  as  error  proxy. 

Key  words'.  Ocean  Ensembles,  Forecasting,  Ocean  Models,  Forecasting  Errors,  Environmental 
Assessment  (43-45N,  9-10. 5E) 


1.  Introduction 

When  considering  numerical  prediction  of  ocean  dynamic  states  using  nested  domains,  several 
sources  of  error  can  contribute  to  cascading  uncertainty  into  state  variable  estimation  (Coelho  and 
Rixen,  2008).  These  sources  of  error  include  the  errors  of  the  initial  and  lateral  boundary 
conditions,  local  forcing,  bathymetry  errors,  numerical  approximations  and  filtering,  errors  due  to 
approximations  when  assimilating  observations,  errors  in  the  forcing  terms  and  un-resolved  scales 
(sub-grid  variability).  To  address  this  problem,  local  unbiased  (correlation)  and  persistent  errors 
(bias)  of  the  Navy  Coastal  Ocean  Modeling  (NCOM)  System  nested  in  global  ocean  domains,  are 
typically  reduced  and  monitored  by  assimilating  dynamical  balanced  analysis  fields  of  state 
variables,  derived  from  observation  networks,  using  the  Navy  Coupled  Ocean  Data  Assimilation 
(NCODA)  system  (e.g.,  Cummings  2005).  This  system  also  provides  an  error  estimate  of  these 
analysis  fields  at  an  analysis  time. 
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In  recent  implementations  (Coelho  et  al.,  2008;  Fabre  et  al.,  2008),  ensemble  based  stochastic 
methods  have  been  used  to  track  these  NCOM  analysis  multi-scale  ocean  errors  by  running  the 
model  several  times  using  different  forcing  and  starting  from  different  initial  conditions.  The 
resultant  ensemble  spread  was  constrained  at  each  new  analysis  time  by  the  new  estimate  of  the 
analysis  errors  using  a  technique  named  Ensemble  Transform  (ET)  (Bishop  and  Tott,  1999).  In 
order  to  be  accurate,  the  perturbed  ensemble  members  should  be  taken  from  a  fairly  large  number 
of  independent  runs  to  resolve  state  variables  error  covariances  and  should  include  all  significant 
sources  of  error  and  uncertainty  (Judd,  et  al.,  2007,  Lermusiaux,  et  al.,  2006).  Since  this  is  not 
easy  to  obtain  in  operational  timeframes,  and  once  a  smaller  number  of  runs  are  selected,  one  can 
expect  the  ensemble  to  perform  differently  inside  the  simulation  domain  and  through  time 
depending  on  the  number  of  the  dominant  error  modes.  This  limitation  motivates  on-going  work 
in  developing  dedicated  metrics  to  diagnose  and  prognoses  ensemble  performances  through  the 
overall  domains  and  forecasting  lead  times. 

In  any  case,  it  is  anticipated  that  a  small  number  of  runs  may  still  provide  useful  information 
under  certain  conditions  (e.g.  when  there  are  no  strong  non-linearity  and  bias  errors  are  on  the 
same  order  of  magnitude  of  the  correlations  errors).  Furthermore,  if  the  ensemble  estimates  define 
a  domain  that  contain  the  most  relevant  features  and  scales  of  the  physical  system,  then  they  can 
be  improved  in  their  consistency  through  calibration  and  post-processing  by  adjusting  their  spread 
and  bias  to  some  training  sequence.  These  methods  have  been  successfully  used  for 
meteorological  ensemble  calibration  (e.g.Doblas-Reyes,  2005;  Hammil,  2007)  and  for  multi¬ 
model  ocean  ensembles  applications  (e.g.  Rixen  and  Coelho,  2007;  Coelho,  2008). 

It  should  be  noted  that  with  a  small  number  of  independent  runs  we  should  not  expect  to  resolve 
the  full  ocean  state  covariances  with  the  original  model  grid  resolution,  but  one  can  expect  a 
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small  number  of  runs  between  10-15  may  still  be  adequate  to  track  single  variable  forecast  errors 
on  a  re-sampled  spatial  domain  as  long  as  the  number  of  independent  variables  can  be  kept  within 
the  order  of  O(103),  following  the  estimates  of  Judd  (2007).  This  note  will  discuss  the  limitations 
of  a  small  ensemble  size  used  during  the  MREA07  trial  and  proposes  a  method  to  improve 
forecast  error  prediction  consistency  for  specific  target  variables,  applicable  also  for  non-state 
variables  estimates  when  there  are  not  many  observations  or  prior  to  use  observations  into  the 
assimilation  process. 

Several  methods  have  been  used  to  perturb  the  initial  conditions  fields  based  on  the  observed 
errors.  In  particular  Bishop  and  Toth  (1999)  proposed  a  technique  named  Ensemble  Transform 
that  allows  computing  dynamically  balanced  initial  conditions  perturbations  that  are  consistent 
with  a  best  estimate  of  the  error  covariance.  On  the  other  hand,  ensemble  calibration  can  also  be 
sought  through  post-processing  using  Bayesian  methods  (e.g.  Gneiting,  2004,  Coelho,  2005  and 
Rixen  and  Coelho,  2005),  within  the  limits  of  the  known  cross-correlations  among  the  observed 
and  modeled  variables.  This  work  combine  both  techniques  as  a  post-processing  method,  applied 
to  local  single  variable  ensemble  spread  calibration.  The  methodology  uses  the  perturbed  model 
statistics  re-scaled  through  an  estimate  of  the  error  variance,  to  obtain  short  term  estimates  of 
posterior  normal  probability  distributions  envelopes  of  a  selected  ensemble  variable. 

The  MREA07  (BP07  and  LASIE  trials),  took  place  off  La  Spezia,  Italy  in  the  spring  and  summer 
of  2007  (e.g.,  LeGac  and  Hermand,  2007).  During  the  trial,  mesoscale  relocatable  NCOM 
implementations  using  the  RELO  system  were  made  available  in  real-time  without  performing 
local  data  assimilation,  though  remote  sensing  and  global  data  was  assimilated  on  the  outer  nests 
used  for  boundary  conditions  and  initialization.  In  standard  implementations  the  RELO  system 
runs  together  with  the  Navy  Coupled  Ocean  Data  Assimilation  (NCODA)  system  that  performs 
observations  quality  control  and  produce  local  analysis  for  assimilation  that  in  the  present  version 
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are  based  on  a  Multi-Variate  Optimum  Interpolation  technique  (e.g.  Cummings,  2005).  NCODA 
also  provides  the  analysis  error  fields  that  are  used  to  re-set  the  ensemble  spread  of  the  initial 
fields  in  operational  ensemble  runs  using  the  same  ET  technique  (e.g.  Fabre  et  al.,  2008).  This 
present  solution  does  not  provide  reliable  analysis  error  covariances  but  it  is  planned  the  NCODA 
system  will  evolve  in  the  near  future  into  using  hybrid  Monte-Carlo  ensembles  (e.g.  Lermusiaux 
et  al.,  2006)  and  Variational  analysis  (e.g.  Nogodock,  et  al.,  2007).  This  will  improve  error 
covariance  estimates  and  produce  analysis  fields  consistent  with  the  boundary  conditions  and 
other  forcing  fields.  For  this  specific  implementation,  the  NCODA  assisted  assimilation  process 
in  the  inner  nests  was  turned  off  to  allow  a  fully  independent  analysis  of  the  model  results  and 
observations,  simulating  a  scenario  where  no  local  data  would  be  available  in  useful  timeframes). 

During  this  tnal  the  free-run  error  fields  of  the  RELO  system  were  estimated  using  an  ensemble 
of  10  independent  runs  with  independent  initial  conditions  starting  from  a  common  field  far  back 
in  time  and  perturbed  through  atmospheric  forcing  using  space-time  deformation  of  the  surface 
forcing  fields  (Xiandong,  2007).  The  ensemble  spread  of  the  free  runs  was  then  re-scaled  in  post¬ 
processing  through  an  Ensemble  Transform  (Bishop  and  Toth,  1999)  using  the  temporal 
variability  as  an  error  proxy.  These  preliminary  error  estimates  were  then  used  for  model 
benchmarking  and  aiming  specific  ocean-acoustic  applications  (e.g.  Carriere  et.al,  this  volume) 
and  to  estimate  the  relative  impact  of  different  observational  strategies  (Coelho  et  al.,  2007). 

2.  RELO-NCOM  Setup 

The  Relocatable  Navy  Coastal  Ocean  Model  (RELO-NCOM)  is  a  scalable,  portable,  and  user- 
friendly  system  for  nowcasting  and  short-term  (2-3  day)  forecasting  simulations  (Rowley,  2007). 
There  are  two  major  components:  1)  NCOM  (Martin,  2000)  and  2)  the  Navy  Coupled  Ocean  Data 
Assimilation  (NCODA)  (Cummings,  2005)  for  data  analysis  and  model  initialization.  For  a  rapid 
configuration,  the  system  relies  on  a  set  of  data  and  products  available  on  a  global  scale 
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(bathymetry,  winds,  analysis  of  the  remote  sensing  data).  These  products  are  generally  on  a  low 
resolution  and  it  is  possible  to  substitute  them  with  local  and  high-resolution  databases.  RELO- 
NCOM  meets  the  naval  requirements  to  generate  real-time  description  of  the  environmental 
variables  and  it  is  operational  at  the  US  Naval  Oceanographic  Office  (NAVO). 

There  is  a  fundamental  difference  between  assessing  an  ocean  model  configuration  in  a  research 
and  an  operational  mode.  Both  need  to  be  designed,  calibrated,  and  evaluated  to  encompass  the 
dominant  dynamics  of  a  given  region.  The  goal  is  to  provide  the  best  possible  representation  of 
the  dynamical  features  of  a  specific  area.  However,  a  predictive  system  that  supports  operational 
applications  must  be  rapidly  relocatable  anywhere  in  the  ocean  (oil-spill  response  and  naval 
operations  are  the  most  relevant  applications),  and  easily  reconfigured.  The  principal  goal  is  to 
provide  good  representations  everywhere  with  the  available  data  (i.e.,  in  spite  of  the  absence  of 
complete  sets  of  observations),  motivating  the  need  to  associate  with  the  system  a  reliable  error 
diagnostics  and  prediction  tool,  to  allow  tracking  consistently  the  error  dynamics. 

For  the  MREA07  trial  the  RELO-NCOM  was  deliberately  set  on  its  default  mode  as  for  a  generic 
application  with  little  or  no  tuning  of  the  physical  and  numerical  parameters.  Furthermore,  no 
MREA07  or  other  data  were  assimilated  into  the  inner  nests.  The  goal  of  this  implementation  was 
to  test  the  modeling  skills  of  these  free  runs  and  estimate  the  relevance  of  the  atmospheric  forcing 
as  a  single  source  of  error. 

The  daily  predictive  cycle  during  MREA07  is  described  as  follows: 

•  NCOM  is  started  from  the  previous  day’s  nowcast  (-24  hr)  and  forced  by  the  available 
operational  winds.  Open  Boundary  Conditions  (OBC)  are  extracted  from  the  simulation 
of  the  parent  domain.  The  OBC  for  the  outer  most  nest  are  extracted  from  NCOM 
configured  on  a  global  scale  at  a  1/8°  resolution  (NCOM-GL)  which  is  operational  the 
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(http://www7320.nrlssc.navy.mil/globalncom/index.html)  (Barron  et  al.,  2006). 
However,  this  procedure  is  not  restricted  to  NCOM-NCOM  nesting;  any  nest  could  be 
coupled  with  several  other  dynamical  model  formulations. 

•  During  the  nowcast,  temperature  (T)  and  salinity  (S)  fields  are  nudged  to  the  nowcast 
fields  of  the  parent  simulations.  The  nudging  during  the  hindcast  phase  has  been 
suggested  to  provide  a  minimum  connection  with  realtime  data  since  NCOM-GL 
assimilates  sea  surface  temperatures  (SST)  and  Modular  Ocean  Data  Assimilation  System 
(MODAS)  synthetics  (with  the  surface  height  derived  from  the  Naval  Layer  Ocean 
Model  (NLOM)  (http://www7320.nrlssc.navy.mil/global_nlom/).  No  data  are  nudged 
after  the  nowcast  (0  hr). 

•  A  short-term  (2-day)  forecast  is  provided.  The  48-hour  interval  has  been  chosen  because 
this  is  the  typical  period  in  which  meteorological  mesoscale  forecasts  are  available  and 
reliable. 

•  The  nested  domains  run  then  in  sequence  using  boundary  conditions  from  the  outer  nests 
(i.e.,  one  way  nesting).  Although  NCOM  provides  a  tile  nesting  approach,  the  default 
procedure  allows  an  easy  and  rapid  configuration  and  assessment  of  each  domain,  and 
more  importantly,  a  possible  different  choice  of  the  vertical  coordinate  between  nests. 
Fig.  1.  illustrates  the  triple  nested  configuration  for  the  MREA07  exercise. 

In  this  model  configuration,  all  domains  are  forced  with  the  Coupled  Ocean  Atmosphere 
Mesoscale  Prediction  System  (COAMPS* +  )  Europe-2  winds  (27km)  (Hodur,  1997)  and  heat 
fluxes  from  0.5°  Navy  Operational  Global  Atmospheric  Prediction  System  (NOGAPS,  Rosmond 
et  al.,  2002).  Monthly  river  discharges  are  extracted  from  the  global  river  data  set  of  NCOM-GL 
(Barron  and  Smedstad,  2002),  with  the  Amo,  Magra,  and  Serchio  transports  provided  by  the 
+  COAMPS  is  a  registered  trademark  of  the  Naval  Research  Laboratory 
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Istituto  Idrografico  Italiano.  The  vertical  resolution  of  each  domain  has  38  a-  and  7  z-levels  (45 
levels).  The  outer  nest  (nestO)  is  at  4km  horizontal  resolution  with  the  primary  purpose  of  serving 
as  a  buffer  zone  between  NCOM-GL’ s  NOGAPS  forcing  and  the  higher  resolution  wind  data  set. 
Nest  I  (2km  resolution)  include  tides.  Tides  are  specified  at  the  boundaries  from  the  Oregon  State 
University  tide  model  (Egbert  and  Erofeeva,  2002).  Nest2  and  nest3  are  at  about  0.6km  resolution 
and  configured  for  the  BP07  (Elba)  and  LASIE07  (LaSpezia)  domains,  respectively.  An  ensemble 
of  10  independent  runs  of  the  inner  nests  was  also  made  available  in  realtime,  using  similar  set¬ 
ups  but  with  perturbed  atmospheric  forcing  using  the  space-time  deformations  method 
(Xiandong,  et  al.,  2007). 

[FIGURE  I] 

One  of  the  most  pressing  issues  of  realtime  operational  forecasting  is  to  provide  the  information 
in  a  timely  manner.  Ocean  forecasts  are  usually  one  of  the  final  components  of  a  long  string  of 
products  developed  at  several  different  centers:  a  delay  in  acquiring  one  of  the  input  data  (e.g., 
winds,  boundary  conditions),  the  classic  computer  breakdowns  (just  to  mention  a  few  issues)  may 
create  a  domino  effect  and  ultimately  a  late  delivery  of  the  forecast.  In  order  to  avoid  delays  in 
the  queue  submission  which  are  often  occurring  at  the  supercomputer  sites,  the  full  forecast  cycle 
is  performed  at  the  Naval  Research  Laboratory  -  Stennis  Space  Centre  (NRLSSC)  on  dual 
processor  Opteron-based  LINUX  platforms.  The  latest  NOGAPS  and  COAMPS  analysis  and 
forecasts  are  usually  available  at  NRLSSC  before  1000GMT,  but  NCOM-GL  daily  hindcasts  and 
forecasts  arrive  at  about  1 130GMT.  Therefore,  to  speed  up  the  delivery  of  the  results,  the  OBC 
for  nestO  are  extracted  from  the  NCOM-GL  72hr  forecast  of  the  previous  day.  This  makes  it 
possible  to  start  the  simulations  at  about  1000GMT  and  complete  the  forecast  cycle  usually 
before  NCOM-GL  latest  files  are  available  at  NRLSSC.  Unfortunately,  only  a  partial  COAMPS 
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data  set  is  archived  at  NRLSSC,  so  the  price  for  this  procedure  is  the  use  of  NOGAPS-O.5  heat 
fluxes. 

The  model  results  are  written  to  NetCDF  files  at  user  specified  z-levels  and  time  increments.  It  is 
important  that  the  z-levels  be  consistent  with  the  NCOM  vertical  grid.  A  coarse  vertical 
resolution  in  the  NetCDF  files  may  remove  features  reproduced  by  the  model;  a  too  fine  vertical 
resolution  increases  the  computational  cost  and  memory  requirement  without  increasing  the 
physical  accuracy  of  the  solutions.  For  this  real-time  exercise,  the  NCOM  fields  were  provided 
on  47-levels  and  at  a  l  hr  increment.  To  reduce  the  amount  of  transferred  data,  only  the  48  hr 
forecast  (i.e.,  no  hindcast)  of  the  model  and  only  a  few  upper  vertical  levels  for  the  ensemble 
spread  were  posted  on  the  MREA07  ftp  server,  generally  at  about  1230GMT  and  1500GMT, 
respectively. 

3.  RELO-NCOM  Control  Analysis  and  data  comparison 

This  note  will  focus  the  analysis  and  discussion  for  the  period  June  10  to  25,  2007  and  for  the  nest 
3  area  only.  In  this  region,  dynamics  were  mostly  dominated  by  a  persistent  cyclonic  gyre 
centered  roughly  at  43  40N  and  9  20W,  modulated  by  smaller  re-circulation  cells  north  and  east, 
closer  to  the  coast.  The  shapes  and  temperature  distributions  of  these  smaller  cells  was  strongly 
perturbed  by  the  wind  forcing.  During  the  “sirocco”  south-easterly  winds  (e.g.  06/19  06:00 
snapshot  displayed  in  Figure  2,  left  panel)  the  average  surface  temperatures  were  higher,  with 
warmer  waters  trapped  closer  to  the  eastern  coast.  During  the  “libeccio”  south-westerly  winds 
((e.g.  06/23  12:00  snapshot  displayed  in  Figure  2,  right  panel),  the  cold  eddy  signature  becomes 
more  noticeable  and  different  recirculation  patterns  can  be  found  between  the  eddy  and  the 
coastline. 

[FIGURE  2] 
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The  Sea  Surface  Temperature  (SST)  images  obtained  from  NOAA  AVHRR  displayed  in  Figure 
3,  although  with  different  resolutions,  concur  with  the  analysis  of  the  previous  paragraph. 

The  water  column  was  strongly  stratified  during  the  whole  period.  Model  temperature  hindcast 
and  forecast  estimates  were  compared  with  160  CTD  profiles  collected  during  the  trial  in  the 
period  June  4-26,  2007  by  three  ships  in  the  area  (RV  Planet,  RV  Leonardo  and  NI  Galatea).  The 
daily  CTOs’  covered  both  deep  and  shallow  water  throughout  most  of  the  surveying  time.  For  this 
work  only  profiles  inside  the  nest  3  domain  were  used.  For  each  CTD,  the  nearest  (in  space  and 
time)  hourly  model  profile  was  extracted.  No  horizontal  or  temporal  interpolation  is  performed  on 
the  model  or  data.  Since  observations  are  on  a  higher  vertical  resolution  relative  to  model 
estimates,  the  model  temperature  at  a  specific  z-level  should  be  compared  with  the  mean  value  of 
the  observed  values  between  the  intermediate  levels  up  and  below  (i.e.  for  the  model  estimate  Tj 
at  level  Zj,  observations  should  be  averaged  between  the  levels  (Zi.j+Zj)/2  and  (Zi+Zj+i)/2).  The 
model  data  comparisons  displayed  in  Figure  4  show  that  temperature  errors  were  more  noticeable 
on  average  at  the  bottom  of  the  well  mixed  layer  (at  roughly  50m  depth),  with  the  surface  waters 
typically  cooler  than  observations  and  warmer  waters  below.  Temperature  errors  were  very  small 
below  the  200m  depth.  It  is  also  noticeable  these  error  characteristics  did  not  change  significantly 
during  the  analysis  period,  though  significant  changes  occur  in  the  forcing  and  dynamic  responses 
as  mentioned  above. 

[FIGURE  3] 

From  these  comparisons  one  can  assume  the  prediction  skills  of  the  model  were  limited,  not 
significantly  above  model  persistency,  such  that  these  free-run  RELO-NCOM  fields  could  be 
considered  as  an  analysis  tool  capable  of  providing  reasonable  spatial  distributions  of  the 
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temperature  fields,  up  to  at  least  48  hours.  This  is  mostly  due  to  the  persistent  nature  of  the 
dominant  local  dynamics  that  did  not  change  significantly  during  the  analysis  period.  In  other 
more  dynamic  areas  one  could  expect  these  free-run  errors  to  increase  significantly  after  a  few 
hours  and  differences  between  forecast  lead  times  also  to  become  more  noticeable. 

Since  there  were  no  significant  differences  between  these  errors,  the  discussion  below  regarding 
error  prediction  will  use  the  0  24  hours  and  24  48  hours  temperature  forecasts  as  equivalent 
estimates. 

4.  Ensemble  re-scaling  using  the  Ensemble  Transform 

The  ocean  is  driven  by  surface  fluxes  that  are  determined  by  the  atmospheric  state  and  are  one 
major  source  of  uncertainty.  Predicted  atmospheric  fields  often  contain  the  forecast  feature  of 
interest,  but  they  can  be  misplaced  in  space  and  time  (e.g.  Hoffman  1995).  This  characteristics 
motivated  attempts  to  represent  forecast  errors  in  terms  of  a  shift  of  a  forecast  in  space  and  time 
similar  to  the  pseudo-random  fields  method  described  by  Evensen  (2003)  and  applied  in  ocean 
ensemble  generation  problems  (e.g.  Demirov,  et  al.,  2003).  For  the  present  work,  the  atmospheric 
forcing  perturbations  used  to  force  the  ocean  ensemble  members  were  produced  using  the  method 
developed  by  Xiandong  et  al.  (2007).  It  uses  only  time  shifts  of  the  forecast,  with  a  choice  of 
parameters  to  provide  a  good  precision  in  the  atmospheric  perturbations,  though  accuracy  may 
not  be  guarantee  over  the  whole  simulation  period. 

[FIGURE  4] 

If  we  neglect  bathymetry,  error  induced  by  numerical  approximations  and  other  sources  of 
possible  model  bias,  the  ensemble  transform  (ET)  method  of  generating  initial  perturbations 
applied  in  atmospheric  ensemble  forecasts  (Bishop  and  Toth,  1999)  can  be  used  to  re-balance  and 


re-shape  the  IC  fields  of  the  ensemble  subset.  Besides  assuring  all  detected  error  growing  modes 
will  be  equally  represented,  the  advantage  of  this  technique  is  such  that:  it  respects  hydrodynamic 
balances  by  ensuring  that  initial  perturbations  are  a  linear  sum  of  forecast  perturbations  from  the 
preceding  forecast;  and  ensures  that  the  initial  perturbations  are  equally  likely  and  orthogonal 
under  a  measure  of  the  probability  of  initial  condition  error  based  on  the  best  available  estimate  of 
initial  condition  error  variance.  This  technique  does  not  provide  though  an  initial  set  of 
background  perturbations  that  need  to  be  introduced  using  complementary  methods,  such  as 
forcing  from  an  ensemble  of  atmospheric  forecasts  as  mentioned  in  the  previous  paragraph. 

As  detailed  in  Toth  and  Bishop  (1999),  through  the  ET  ensemble  generation  technique,  K  forecast 
perturbations  of  N  state  variables  X°  ( NxK ),  can  be  transformed  into  a  set  of  perturbations  Xr 
that  are  consistent  with  the  background  error  analysis  covariance  P"  ,  using 

Xr  =  X°T 

where  T  is  a  transformation  matrix  determined  by  the  eigenvectors  and  eigenvalues  of  the 
projections  of  the  magnitude  of  the  predicted  analysis  perturbations  on  the  inverse  of  the  error 
analysis  covariance  matrix.  If  the  number  of  ensemble  members  equals  the  number  of  state 
variables,  this  projection  guarantees  the  perturbations  covariance  to  be  equal  to  the  error 
covariance. 

Through  this  transform  we  can  then  obtain  a  set  of  perturbed  fields  that  are  consistent  with  an 
independent  estimate  of  the  error  covariance.  In  operational  implementations  these  initial  fields 
are  used  as  new  initial  conditions  for  the  K  independent  ensemble  runs,  providing  a  method  to 
assimilate  the  observed  errors  into  the  ensemble  forecasts.  For  the  present  application  and  to  use 
this  method  in  post-processing  a  persistency  assumption  during  the  48hours  forecast  cycles  is 
taken,  regarding  the  projection  of  the  ensemble  covariances  into  the  observed  errors. 
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5.  MREA07  Error  Predictions 


For  the  present  application  since  no  data  is  to  be  used  the  ET  is  computed  using  the  temperature 
48  hours  forecast  time  variances,  as  estimated  by  the  RELO-NCOM  free  runs,  producing  a 

diagonal  error  covariance  matrix  P* .  Besides  allowing  for  a  faster  transform,  this  approach 

allows  to  keep  the  shapes  of  the  off-diagonal  terms  (spatial  cross-correlations)  as  estimated  by  the 
ensemble,  while  consistently  re-scaling  the  analysis  errors,  without  introducing  further  analytical 
or  numerical  approximations. 

The  temperature  estimates  ensemble  spatial  correlations  are  then  updated  only  by  the  RELO- 
NCOM  independent  runs.  This  method  allows  keeping  error  covariance  updates,  without  the  cost 
of  computing  and  inverting  very  large  matrices.  Furthermore,  since  only  a  limited  number  of 
ensemble  members  are  available,  this  method  limits  the  growth  of  spurious  cross-correlations. 
The  same  transform  matrix  T  is  applied  to  all  time  steps  of  the  ensemble  estimates. 

The  resulting  ensemble  spread  (standard  deviations)  for  each  temperature  estimate  are  then 
compared  against  the  absolute  value  of  the  RELO-NCOM  vs.  data  mismatches  and  displayed  in 
scatter  diagrams  as  those  shown  in  Figure  5  for  days  Jun  13  and  14,  before  and  after  applying  the 
ET.  The  statistical  significance  of  each  of  these  individual  estimates  (small  blue  dots)  is 
negligible,  such  that  they  are  grouped  in  equally  populated  bins  with  1000  elements,  defined 
along  the  ensemble  spread  axis.  These  bins  displayed  inside  the  scatter  diagrams  as  large  red  dots 
will  have  similar  likelihoods  and  will  be  statistically  relevant.  For  the  ensemble  to  be  accurate, 
bins  should  be  aligned  along  the  main  diagonal,  highlighted  as  a  black  line  on  the  plots.  The 
green  rectangles  around  the  bins  show  the  standard  deviations  of  each  bin  along  each  axis  (error 
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and  ensemble  spread).  Other  relevant  statistic  is  the  mean  ratio  between  measured  error  vs. 
ensemble  spread,  (Err/Std  in  the  figures)  that  should  be  close  to  1  for  the  ensemble  to  be  accurate. 

[FIGURE  5] 

The  graphics  in  Figure  5  left  of  the  black  line  show  the  scatter  diagrams  for  days  13  (left  upper 
plot)  and  day  14  (right  upper  plot)  computed  from  the  ensemble  before  post-processing.  From  the 
bin  distribution  we  can  see  the  ensemble  to  have  a  positive  spread-skill  relationship,  through  all 
ranges  of  the  observed  errors,  such  that  estimates  of  smaller  ensemble  spread  are  well  correlated 
with  smaller  errors  and  estimates  of  larger  error  are  well  correlated  with  the  larger  errors,  through 
all  ranges  of  observed  errors.  However,  we  can  see  the  ensemble  was  grossly  under-predicting  the 
magnitudes  of  the  observed  errors  in  roughly  one  order  of  magnitude.  This  is  most  likely  due  to 
the  fact  the  initial  fields  and  other  major  sources  of  error  besides  atmospheric  forcing  were  not 
being  properly  perturbed. 

The  data  of  June  13  was  used  as  the  initial  day  to  start  the  procedure  and  adjust  the  ensemble 
spread  to  the  observed  error.  For  this  purpose,  a  multiplication  factor  of  4  was  estimated  from  the 
data  and  applied  to  the  temporal  standard  deviations  used  to  compute  the  ET  throughout  the 
simulation  period.  This  value  was  estimated  iteratively  in  order  to  bring  the  ratio  Err/Std  from  a 
value  of  1 1  before  the  transform  to  1.  As  a  result,  the  red  bins  also  became  closer  to  the  main 
diagonal  as  we  can  see  on  the  scatter  diagrams  right  of  the  vertical  black  line  in  Figure  5.  For  the 
following  day  represented  by  the  24-48  hour  forecast  this  ratio  increased  slightly  to  1.5,  though 
the  bins  remained  close  to  the  main  diagonal. 

Other  relevant  result  from  Figure  5  is  the  spatial  distribution  of  the  error  estimates.  In  the  lower 
color  maps  one  can  see  the  ensemble  spread  at  the  surface  for  days  13  and  14.  The  black  crosses 
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show  the  points  were  data  was  collected  during  those  days  respectively.  One  can  see  the  spatial 
patterns  were  not  strongly  changed  by  the  transform  and  the  areas  with  larger  estimated  errors  are 
shaped  along  the  boundaries  of  the  persistent  cyclonic  eddy  in  the  SW  portion  of  the  domain  as 
one  could  expect.  The  sampling  locations  during  these  two  days  included  several  runs  across  the 
boundaries  of  this  cyclonic  gyre. 

Since  the  ET  was  using  the  temporal  standard  deviation  to  re -scale  the  ensemble  spread  one  could 
argue  that  the  information  contained  in  the  ensemble  would  be  erased  and  time  variability  would 
be  the  dominant  error-proxy.  In  order  to  evaluate  this  hypothesis  the  same  scatter  diagrams  were 
computed  using  the  temporal  standard  deviation  instead  of  ensemble  spread,  as  displayed  in 
Figure  6.  To  keep  an  equivalent  accuracy  a  multiplication  factor  of  7.8  was  also  applied  to  set  the 
ratio  Err/Std  to  1  for  the  day  13  data.  From  the  scatter  diagrams  one  can  see  this  error  proxy  keeps 
similar  positive  spread-skill  relations,  though  the  spatial  distribution  of  errors  is  significantly 
different  from  those  estimates  by  the  ensemble  and  not  so  well  correlated  with  the  dominant 
dynamics. 

Using  the  tuning  parameters  estimated  for  day  13,  one  can  estimate  the  ensemble  spread  and  the 
time  variability  error  proxy  for  the  following  forecast  days.  Since  observations  were  made  until 
June  25,  Figure  7  display  the  same  diagrams  for  the  last  two  days  of  June  24  (0/24  hours  in  the 
labels)  and  25  (24/48  hours  in  the  labels)  when  model-data  comparisons  were  possible.  The  four 
plots  panel  in  the  left  shows  the  results  using  the  transformed  ensemble  and  the  panel  in  the  right 
shows  the  same  results  using  the  time  variability  proxy.  One  can  see  the  ensemble  spread  was 
kept  consistent  with  the  dynamics  and  the  performance  of  both  the  transformed  ensemble  and 
time  variability  as  error  proxy  seem  close  in  performance.  However,  looking  to  the  spatial 
distribution  of  the  predicted  surface  temperature  errors  as  displayed  in  the  lower  color  maps  for 
days  June  24  and  25  one  can  see  the  ensemble  responded  consistently  with  the  “Sirocco’’  and 
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“Libeccio”  wind  events,  spreading  the  areas  of  larger  uncertainty  around  the  cyclonic  eddy,  not  so 
well  represented  by  the  time  variability  proxy. 

[FIGURE  6] 

In  order  to  obtain  more  objective  performance  estimates,  daily  performance  statistics  were 
computed  as  displayed  in  table  1.  These  include  the  ration  Err/Std  as  an  estimate  of  the  error 
estimate  accuracy,  the  bins  correlation  coefficient  (C)  as  an  estimate  of  the  spread-skill  and  the 
bin  deviation  from  the  main  diagonal  (Bin  Bias  -  BB)  as  an  estimate  of  the  error  estimates  bias. 

[FIGURE  7] 

Overall,  during  the  period  June,  13  to  25  the  positive  spread-skill  was  kept  for  all  estimates 
(ensemble  with  and  without  transform  and  time  variability),  with  the  ensemble  performing 
slightly  better  showing  a  0.8  correlation  coefficient  among  the  bins  while  the  time  proxy  had  a  0.7 
coefficient.  The  ratio  Err/Std  was  also  kept  consistently  through  this  period  such  that  on  average 
through  this  period  the  ensemble  value  was  13.4,  the  ET  was  kept  as  1  and  the  time  proxy  as  1.1. 

The  mean  differences  between  bin  coordinates  (i.e.  deviations  from  the  main  diagonal)  can  also 
be  used  as  an  error  bias  estimate.  Through  this  12  days  period  (June  13  to  25)  the  ensemble 
estimates  after  the  transform  remained  unbiased  while  the  original  ensemble  had  a  value  of  0.4 
and  the  time  variability  proxy  showed  also  a  negligible  negative  bias  of  0.03. 

6.  Concluding  Remarks 

The  work  presented  above  showed  that  some  level  of  predictability  of  stochastic  environmental 
variables  through  numerical  modeling  could  be  achieved  using  Monte-Carlo  methods,  producing 
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ensemble  based  error  estimates  along  with  the  predicted  state  variables,  even  using  a  limited 
number  of  ensemble  runs.  However,  the  system  performance  will  be  space  and  time  dependent 
requiring  an  accurate  metrics  system  to  produce  both  diagnostics  and  prognostics  of  the  precision 
and  accuracy  of  the  outputs. 

The  Ensemble  Transform  (ET)  approach  was  successfully  applied  for  free-run  ocean  Mesoscale 
error  prediction  calibration,  by  re-scaling  RELO-NCOM  ensembles  produced  through 
atmospheric  perturbations.  Independent  data  was  used  for  this  analysis  where  the  model  runs 
were  not  assimilating  any  local  data.  Results  show  the  ensemble  spread  did  not  diverge  and  was 
consistent  with  the  observed  dynamics  throughout  the  simulation  period.  The  ensemble  showed  a 
positive  spread-skill  through  all  ranges  of  the  observed  errors. 

Comparisons  of  ensemble  spread  of  the  temperature  profiles  with  local  observed  errors  and  time 
variability  (assumed  as  an  error  proxy)  showed  they  were  consistent  through  a  12  days  analysis 
period.  The  ET  calibrated  ensemble  had  slightly  better  performance  statistics  then  the  time- 
variability  error  proxy,  most  likely  due  to  the  fact  the  ensemble  predicted  errors  were  better 
correlated  with  the  local  observed  dynamics. 

Results  show  the  ensemble  spread  did  not  diverge  and  was  consistent  with  the  observed  dynamics 
throughout  the  simulation  period.  Furthermore,  comparisons  of  ensemble  spread  of  the 
temperature  profiles  with  local  observed  errors  and  time  variability  (assumed  as  an  error  proxy) 
showed  they  were  consistent  through  the  12  days  analysis  period,  with  performances  above  the 
non-calibrated  ensemble  estimates  and  time-variability  used  as  error  proxy.  Overall  error 
estimates  became  unbiased  and  the  system  was  able  to  accurately  separate  large  errors  from 
smaller  errors  with  a  positive  spread-skill  relationship,  through  all  ranges  of  the  observed  errors. 
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Table  1  -  This  table  shows  the  daily  mean  values  of  the  ration  between  individual  observed  error 
magnitudes  vs  the  correspondent  ensemble  standard  deviation(Err/Std),  the  correlation  coefficient 
or  linear  regression  slope  of  the  1000  point  bin  averages  (Corr.Coef)  and  the  difference  between 
the  bins  ensemble  standard  deviation  and  bin  errors  in  degrees  C  (Bin  Bias  BB).  Each  one  of 
these  estimates  was  computed  for  the  ensemble  withouth  post-processing  (Ens),  with  the  ET  post¬ 
processing  (ET)  and  for  the  post-processed  time  variability  used  as  an  error  proxy  (Time).  The 
raw  at  the  bottom  shows  the  overall  averages  during  the  experiment. 


Day 

Err/Std 

Corr.Coef. 

Bin  BIAS  (BB) 

Ens 

ET 

Time 

Ens 

ET 

Time 

Ens 

ET 

Time 

06/13 

11.0 

1.0 

1.0 

0.84 

0.84 

0.75 

0.3 

0.0 

0.0 

06/14 

16.3 

1.5 

1.5 

0.74 

0.73 

0.67 

0.4 

0.1 

0  1 

06/17 

20.3 

1.6 

1.8 

0.79 

0.80 

0.55 

0.3 

0.1 

0.2 

06/18 

10.0 

0.8 

1.0 

0.67 

0.67 

0.89 

0.5 

-0.2 

0.0 

06/20 

17.4 

1.4 

1.8 

0.48 

0.53 

0.38 

0.4 

0.1 

0.2 

06/21 

12.1 

0.9 

1.0 

0.89 

0.90 

0.85 

0.3 

0.0 

0.0 

06/22 

10.3 

0.8 

0.9 

0.85 

0.86 

0.76 

0.3 

-0  1 

-0.1 

06/23 

13.0 

1.0 

0.9 

0.90 

0.91 

0.90 

0.3 

0.0 

0.0 

06/24 

11.9 

0.8 

0.6 

0.85 

0.85 

0.94 

0.2 

-0.1 

-0.2 

06/25 

11.9 

0.8 

0.7 

0.70 

0.69 

0.46 

1.3 

-0.3 

-0.5 

MEAN 

13.4 

1.0 

1.1 

0.8 

0.8 

0.7 

0.43 

0.00 

-0.03 
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Fig.l.  The  triple  nest  configuration  for  MREA07. 

Fig.  2  -  RELO-NCOM  upper  layer  temperature  snapshots  for  the  days  06/19  (left  panel)  and 
06/23  (right  panel).  The  shapshots  hours,  displayed  in  the  images,  correspond  to  the  wind 
maximum  stress  for  each  day.  During  the  19th  winds  were  predominantly  south-easterly 
(“Sirocco”)  and  during  the  23rd  they  were  predominantly  south-westerly  (“Libeccio”).  Both 
panels  display  how  flow  patterns  changes  around  the  persistent  gyre  in  the  South-West  comer, 
with  warmer  waters  intruding  northward  during  the  “Libeccio”  event. 


Fig.  3  -  NOAA  AVHRR  Sea  Surface  Temperature  estimates  for  06/19  (left  panel)  and  06/23 
(right  panel).  During  the  19th  winds  were  predominantly  from  the  south-east  (“Sirocco”)  and 
during  the  23rd  from  the  south-west  (“Libeccio”).  Images  were  produced  by  automatic  processing 
using  NURC  TERASCAN  software. 

Fig.  4  -  RELO-NCOM  water  temperature  bias  and  RMS  error  estimates.  The  four  panels  in  the 
left  show  the  RMS  errors  along  each  simulation  day  (24  hours  period),  using  different  model 
estimates  compared  with  the  observations.  The  color  plot  named  “A04”  in  the  upper  left  uses 
hindcast  atmospheric  forcing  fields,  the  plot  named  “Pers”  uses  model  persistency  (hour  0 
snapshot)  and  the  plots  below  named  “F24”  and  “F48”  use  24  and  48  hours  lead  forecasts 
respectively.  The  four  panels  in  the  right  show  the  error  bias  (24  hours  mean  errors)  using  the 
same  model  estimates. 

Fig.  5  Error  scatter  plots  computed  using  the  run  of  June  13.  The  upper  scatter  diagrams  show 
the  ensemble  spread  vs.  observed  forecast  error  before  re -scaling  (Figure  5-a)  and  after  re -scaling 
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(Figure  5-b).  The  forecast  errors  were  computed  using  the  0-24  hour  forecasts  (panels  in  the  left) 
and  using  the  24-48  hour  model  forecasts  (panels  in  the  right).  The  color  plot  below  each  scatter 
diagrams  show  the  surface  temperature  error  estimate  (ensemble  standard  deviation)  at  hour 
00.00  (left)  and  24:00  (right)  relative  to  the  simulation  day  and  the  white  crosses  depicts  the 
locations  used  for  model-data  comparison. 

Fig.  6  -  Same  as  Figures  5-b,  but  using  the  time  variability  as  an  error  proxy  instead  of  the 
ensemble  spread  as  an  error  estimate 


Fig.  7  Same  results  as  described  for  Figure  5-b  (on  the  left)  and  Figure  6  (on  the  right)  but  for 
the  model  run  of  June  24.  The  panels  left  of  the  vertical  line  show  the  results  using  the  calibrated 
ensemble.  Panels  in  the  right  show  the  same  results  but  using  the  time  variability  as  an  error 
proxy. 
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Abstract 

Multi-model  Super-Ensembles  (SE)  which  optimally  combine  different  models,  have 
been  shown  to  significantly  improve  atmospheric  weather  and  climate  predictions. 
In  the  highly  dynamic  coastal  ocean,  the  presence  of  small-scales  processes,  the  lack 
of  real-time  data,  and  the  limited  skill  of  operational  models  at  the  meso-scale  have 
so  far  limited  the  application  of  SE  methods  for  acoustic  Rapid  Environmental 
Assessment  purposes.  In  the  framework  of  the  BP07  experiment  conducted  South 
East  of  Elba,  sound  speed  prediction  skills  of  various  SE  techniques  combining 
operational  model  outputs  and  in-situ  measurements  arc  assessed.  Results  suggest 
that  SE-based  predictions  provide  more  robust  24hr  forecasts.  A  detailed  acoustic 
propagation  sensitivity  study  at  different  frequencies  and  ranges  also  reviews  the 
potential  of  these  predictions  for  acoustic  inversion  and  tomography  efforts. 

Keywords:  Oeean-aeoustie  predictions,  multi-model  super-ensemble,  sound  speed, 
Kalman  filter,  data  assimilation,  tomography 

1.  Introduction 

An  increasing  number  of  models  are  routinely  providing  operational  (atmospheric) 
weather  forecasts  and  climate  predictions  (Palmer,  2004)  but  prediction  skill  is  inherently 
limited  for  a  number  of  reasons,  including  simplifications  in  physical  processes,  errors  in 
initial  conditions  and  boundary  conditions,  numerical  schemes,  etc.  The  use  of  data 
assimilation  techniques  ( e.g .  Bennett,  1992;  Wunseh,  1996;  Robinson  et  al,  2004; 
Bennett,  2002;  Evensen,  2006)  to  regularly  correct  for  model  drifts  may  compensate  to 
some  extent  for  loss  of  predictability  with  time  (Lorenz,  1963).  Model  ensembles  have 
become  an  important  means  of  investigating  dispersion  problems  (Galmarini  ct  al  2001, 
2004),  tracking  individual  model  errors,  increasing  forecast  skill,  and  reducing 
uncertainties  (Lermusiaux  et  al,  2006)  in  highly  dynamic  and  complex  environments 
where  predictability  is  limited.  The  multi-model  Super-Ensemble  (SE)  technique 
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(Krishnamurti  et  al  1999,  2000a,  2000b),  which  uses  an  optimized  combination  of  an 
ensemble  of  models  has  previously  been  demonstrated  to  improve  weather,  seasonal  and 
interannual  forecast  skill  in  atmospheric  (Kumar  et  al  2003;  Shin  and  Krishnamurti 
2003a;  Yun  et  al  2005;  Mutemi  et  al  2007)  and  ocean  (Logutov  and  Robinson  2005; 
Rixen  and  Ferreira-Coelho  2005,  2007;  Rixen  et  al,  2007,  2008)  models  over  simple- 
ensemble  and  bias-removed  ensemble  means.  SE  methods  (Williford  et  al  2003)  have 
been  further  improved  by  the  use  of  dynamic  (Shin  and  Krishnamurti  2003b,  Rixen  et  al 
2008),  regularization  (Yun  et  al  2003),  non-linear  (Rixen  and  Ferreira-Coelho  2007)  and 
probabilistic  (Rajagopalan  et  al  2002)  techniques.  These  methods  all  aim  at  finding  a 
combination  of  models  that  optimally  agrees  with  reference  data  over  a  training  period 
(the  hindeast,  regression  or  fit);  this  combination  is  subsequently  used  to  produce  a  SE 
forecast  obtained  by  weighting  individual  model  forecasts.  A  critical  aspect  for  all  super- 
ensemble  methods  is  therefore  whether  the  regression  solution  is  capable  of  extrapolation 
in  time  and  is  applicable  to  future  events.  In  other  words,  is  the  learning  adequate  to 
provide  generalization  skills? 

Operational  implementation  of  SE  methods  in  Numerical  Weather  Prediction  (NWP) 
centers  is  quite  straightforward  due  to  the  reliability  of  observational  data  streams  and  the 
robustness  of  the  models.  On  the  other  hand,  in  the  ocean,  the  lack  of  long  real-time  data 
time  series  -  especially  in  shallow  waters  -  and  a  limited  suite  of  operational  models 
have  so  far  limited  the  application  of  such  promising  techniques  in  an  operational 
framework.  The  limitations  for  in-situ  observations  in  the  coastal  and  shallow  water 
environment  are  mostly  due  to  heavy  maritime  traffic,  intense  fishing  activity  and 
mechanical  and  biological  stress  on  sensors  and  platforms. 

A  pioneering  study  was  conducted  during  the  MREA04  (Maritime  Rapid  Environmental 
Assessment)  field  experiment  along  the  Portuguese  coasts  to  investigate  the  potential 
benefit  of  SE  techniques  for  acoustic  purposes  and  concluded  that  simple  linear- 
regression  based  multi-model  prediction  were  able  to  improve  significantly  sound  speed 
prediction  skills  at  24hr  lead  time  (Rixen  and  Ferreira-Coelho,  2004). 

Real-time  oeean-aeoustie  predictions  and  data  assimilation  can  be  very  useful  but  require 
a  precise  understanding  of  the  full  transfer  of  uncertainties  from  the  ocean  to  the  acoustic, 
e.g.  using  ensemble  schemes  (Lermusiaux  et  al  2002;  Lermusiaux  and  Chiu,  2002; 
Robinson  and  Lermusiaux,  2004;  Lermusiaux  et  al,  2006.  Adaptive  sampling  schemes 
and  their  impact  on  acoustic  propagation  may  help  reducing  uncertainties  in  specific 
regions  of  interest  (e.g.  Heaney  et  al,  2007;  Wang  et  al,  2008  this  issue;  Yilmaz  et  al, 
2008). 

In  the  framework  of  the  BP07  field  experiment  conducted  South-East  of  Elba  in  Spring 
2007  (Le  Gae  and  Hermand,  2007),  we  investigated  the  use  of  dynamic  SE  techniques 
based  on  the  Kalman  filter  (Kalman,  1960)  to  allow  for  a  temporal  evolution  of  model 
combinations,  which  form  the  basis  of  the  present  work,  described  in  section  2.  A 
thorough  acoustic  propagation  sensitivity  study  is  carried  out  in  section  3  to  assess  the 
potential  of  SE  predictions  for  acoustic  inversion  and  tomography  purposes.  Specifically, 
multi-frequency  correlations  and  uncertainties  are  investigated  in  detail.  In  a  companion 


paper  (Carriere  et  al,  2007,  this  issue),  full-field  tomography  and  Kalman  tracking  of  the 
range-dependent  sound  speed  is  investigated  for  the  same  field  experiment. 


2.  The  BP07  experiment:  ocean  observations  and  predictions 

An  intense  joint  ocean-acoustic  observational  program  and  prediction  effort  took  place 
during  the  BP07  experiment  and  is  described  in  details  in  the  field  trial  report  (LcGac  and 
Hermand,  2007)  and  in  a  number  of  companion  publications  (Carriere  ct  al,  Meyer  ct  al, 
Lam  et  al,  Coclho  ct  al,  this  issue). 


The  main  focus  of  the  experiment  was  on  a  small  area  Southeast  of  Elba  Island,  Italy  and 
in  particular  on  a  section  in  the  middle  of  area  ‘REA  I’  (AB  transect  in  Fig.  1)  for  which  a 
detailed  gcoacoustic  characterization  of  the  seafloor  and  subscafloor  is  available 
following  the  Yellow  Shark  94  inversion  results  (Hermand  and  Gcrstoft  1996,  Hermand 
1999).  However  the  ocean  monitoring  effort  was  not  entirely  focused  on  the  A-B  transect 
but  also  covered  the  wider  BP3  area  in  general. 
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Fig.  1.  MREA/BP07  test  areas.  The  boxes  for  ocean  monitoring  and  prediction 
(BP3)  and  for  acoustic  characterization  (REA1)  are  shown.  Stars  in  yellow,  green 
and  gray  show  CTD  casts  collected  by  R/V  Aretusa  (Italian  Navy),  R/V  Leonardo 


(NURC)  and  R/V  Snellius  (Royal  Netherlands  Navy),  respectively,  during  the  BP07 
time  frame,  April  16  -  May  4,  2007.  Section  AB  has  been  one  of  the  focuses  of  the 
ocean-acoustic  experiments. 

Regular  CTD  were  collected  by  three  vessels  and  provided  a  reasonable  spatio-temporal 
coverage.  In  addition,  two  thermistor  strings  15km  apart  and  equipped  with  1 1  sensors, 
were  deployed  by  NRV  Leonardo  at  position  A  and  B  from  19  April  to  1  May.  They 
eovered  water  depths  respectively  from  13.3m  to  63.5m  and  from  13.5m  to  53.5m  with 
spaeing  of  5m.  Sampling  rate  was  set  to  2  minutes.  Sensors  7  and  1 1  at  mooring  A  failed 
and  were  withdrawn.  Note  (Fig.  2)  the  high  temporal  variability  (some  patterns  have 
cycles  shorter  than  an  hour)  and  the  strong  differences  between  the  two  moorings. 
Surfaee  heating  is  obvious  at  station  B,  whilst  station  B  shows  the  presence  of  a  more 
complex  mixed  layer  on  two  occasions  and  a  more  cyclic  pattern  in  temperature 
evolution. 


Fig.  2  Raw  temperature  thermistor  string  data  time  series  over  depth  (m)  collected 
respectively  at  the  A  and  B  endpoints  of  the  transect. 


These  measurements  were  complemented  by  selected  MVP  surveys  (see  Lam  et  al,  this 
issue),  oceanographic  data  from  acoustic  drifting  buoys  and  profiles  collected  from  R/V 
Snellius  launch  or  rubber  boat  (not  shown  here). 

All  data  collected  in  the  vicinity  of  the  AB  transect  represent  a  pool  of  1 8866  temperature 
data  and  6882  salinity  data,  which  were  objectively  analyzed  and  gridded  using  an 
Optimal  Interpolation  (01)  technique  ( e.g .  Bretherton  et  al,  1976;  Rixen  et  al,  2001)  with 
an  horizontal  grid  resolution  of  715m,  a  vertical  resolution  of  lm  and  a  temporal 
resolution  of  3hrs.  The  3D  box  over  section  AB  for  period  16  April-3  May  2007  was 
spanning  15km,  110m  and  13  days.  To  compensate  for  sensor  noise  and  sensor  inter- 
calibration  issues,  the  noise-to-signal  ratio  was  set  to  1  after  cross-validation  of  this 
parameter.  Because  of  the  amount  of  data,  intractable  for  direct  optimal  interpolation 
techniques,  the  domain  of  analysis  was  split  into  sub-domains  larger  than  typical 
correlation  lengths  to  overcome  the  computational  burden,  a  technique  known  as  sub- 
optimal  interpolation.  The  spatial  correlation  lengths  were  set  to  4m  on  the  vertical  and 
2km  on  the  horizontal.  The  temporal  correlation  length  was  set  to  1  day.  Lower  values 
were  creating  unphysical  results  because  of  the  non-uniform  distribution  of  observations 
in  space  and  time.  A  background  field  was  obtained  by  spatio-temporal  linear  regression 
over  the  whole  experimental  period  (also  known  as  First  Guess  At  Appropriate  Time  - 
FGAT)  and  the  OI  analysis  was  computed  on  resulting  anomalies.  Multivariate  analysis 
techniques  were  not  explored  here.  Sound  speed  fields  were  derived  from  the  T/S 
analyses  (Fofonoff  and  Millard,  1983). 

The  sound-speed  field  over  section  AB  exhibits  both  spatial  and  temporal  variability  as 
illustrated  in  Fig  3  and  4.  Strong  heating  can  be  observed  on  the  upper  10m  at  the  surface 
over  the  period  resulted  in  a  strong  pycnocline  and  high  sound-speed  vertical  gradients  at 
around  15m  depth..  The  area  was  subject  to  strong  mesoscalc  activity  as  well  (see  also 
Carriere  et  al,  this  issue),  as  illustrated  in  the  two  high  sound  speed  surface  cells  (Fig.  4). 
The  diurnal  cycle  has  been  smoothed  out  by  the  analysis  which  was  not  able  to  preserve 
all  the  spectral  information  because  of  the  limited  amount  of  data. 
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Fig  3.  Sound  speed  (m/s)  on  section  AB  versus  depth  for  the  start  (top  left)  to  the  end 
(bottom  right)  of  the  experiment. 
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Fig  4.  Sound  speed  (m/s)  temporal  evolution  over  depth  at  different  latitudes  from 
location  A  (top  left)  to  B  (bottom  right). 

2.2  Models  and  adaptive  sampling 

A  specific  prediction  system  was  set  up  at  NRLSSC  for  the  MREA/BP07  sea  trial.  As 
described  in  details  in  Coelho  et  al  (this  issue),  it  was  built  around  the  NCOM  model. 
Starting  from  the  global  NCOM  model  at  4-km  resolution,  two  nested  models  were  set  up 
in  order  to  cover  the  full  MREA07  area  respectively  with  a  2-km  coarse  resolution 
(NCOM  COARSE)  and  a  0.6  km  fine  resolution  (NCOM  FINE).  SST  data  and  MODAS 
synthetics  were  indirectly  assimilated  in  the  models  through  the  Global-NCOM  models, 
while  the  nested  ones  directly  assimilated  COAMPS-Europe  2  wind  forcing  and 
NOGAPS  heat  fluxes.  48-h  forecasts  were  made  available  on  a  daily  basis. 

NCOM  ensemble  runs  were  used  to  minimize  forecast  error  covariance  by  exploring 
various  observational  patterns  (see  Coelho  et  al,  this  issue).  Error  covariance  analyses 
were  transmitted  by  NRLSSC  to  NURC  and  two  CTD  surveys  were  specifically  designed 
for  this  purpose  and  used  in  a  parallel  for  the  MSEAS-HOPS  ocean  modeling  effort  with 
a  sensitivity  study  of  acoustic  propagation  and  probability  of  detection  at  low  frequencies 
(sec  Lam  et  al  for  details,  this  issue).  The  proposed  CTD  sampling  strategy  was  then 
adapted  according  to  real-time  on-site  constraints  (ship  time  available  in  between 
acoustic  and  geoaeoustie  runs,  weather  and  sea-state,  etc). 


2.3  Super-ensemble  predictions 

The  analysis  field  and  NCOM  model  outputs  were  used  to  explore  various  ensemble 
prediction  techniques. 

The  simple  ensemble  mean  (hereafter  ENSMEAN)  does  not  use  observations  over  the 
training  period  and  thus,  cannot  really  be  considered  as  a  SE  technique.  However  it  is 
also  a  used  method,  since  it  is  usually  expected  to  provide  better  foreeasts  than  individual 
models  (Kalnay  and  Ham,  1 989). 

The  unbiased  ensemble  mean  (hereafter  UNBIASED  ENSMEAN)  corrects  for  biases  on 
each  individual  model,  based  on  observations  during  a  training  period.  These  unbiased 
models  arc  then  averaged. 

The  linear  regression  SE  technique  (LINREG)  consists  in  finding  a  linear  combination  of 
the  models,  minimizing  (in  the  least  mean  squares  sense)  its  departure  from  observations 
during  a  training  period.  The  resulting  weights  are  then  used  to  combine  numerical 
foreeasts.  This  method  can  be  improved  by  normalizing  models  and  adding  a  constant 
model  (i.e.  bias  or  independent  term),  hereafter  appended  with  suffix  NORM). 
Col  linearities  between  the  models  can  also  be  removed  by  retaining  only  a  certain 
percentage  of  variability  of  the  models  by  applying  an  Empirical  Orthogonal  Function 
(EOF)  (also  known  as  Principal  Component  Analysis  -  PCA)  on  the  models  -  in  the 
present  study  95%  of  the  variance,  with  suffix  EOF  -  which  results  in  fewer  models  and 
improves  the  generalization  capabilities  of  the  SE.  This  method  is  hereafter  referred  to  as 
LINREGNORMEOF. 

These  techniques  are  well  known  and  have  been  tested  in  various  oeeanographie  contexts 
(Rixen  and  Ferreira-Coelho  2005,  2007,  Rixen  et  al  2007,  2008;  Lenartz  et  al,  this  issue; 
Vandcnbuleke  et  al,  submitted).  Training  period  is  chosen  a  priori  and  all  observations 
are  equally  important.  Naturally  however,  more  reeent  data  should  be  more  relevant  and 
weights  should  be  adapted  according  to  reeent  model  skills. 

Sequential  data  assimilation  techniques  ean  continuously  adapt  the  weights  during  a 
training  period  when  observations  are  available  up  to  the  present  time  when  weights  are 
frozen  and  used  to  eombine  available  foreeasts  for  the  future.  The  recursive  Kalman  filter 
(Kalman,  1960)  solution  is  well  suited  for  this  purpose  and  is  briefly  described  in  the 
context  of  SE.  It  consists  of  two  eonseeutive  steps. 

1)  prediction  step: 


*'(/,)  =  A/,  *•(/,_,) 


(1) 

(2) 


2)  eorreetion  step: 


K  =  Pf(ti)HT[R  +  HPf(ti)HTy' 
x°{ti)  =  xf{ti)  +  K[y0-Hxf{ti)\ 
Pa{t,)  =  Pf(ti)-KHPf{ti) 


(3) 

(4) 

(5) 


The  state  vector  x  with  covariance  P,  contains  the  weights  on  the  models  in  the  SE 
combination.  Superscript / denotes  forecast  state  after  prediction  steps;  and  superscript  a 
stands  for  analyzed  state  after  the  correction  steps  using  observations.  The  state  vector  x 
is  initialized  with  a  best  guess  obtained  from  the  LINREG  solution.  The  P  is  set  to  0.5 
initially  as  we  expect  P  to  be  far  off  the  optimal  value  at  the  beginning  of  the  training. 

In  the  context  of  SE,  the  model  matrix  M,  with  error  covariance  Q,  is  the  identity  matrix. 
Standard  deviations  of  the  model  error  for  individual  weights  are  set  to  the  variability  of 
weights  for  LINREG  solutions  of  various  short  sub-training  periods,  providing  the 
diagonal  terms  for  the  error  covariance  Q. 

Observations  arc  represented  in  the  vector  y,  with  error  covariance  R.  In  the  present 
study,  y  is  the  analyzed  sound  speed,  obtained  by  optimal  interpolation  of  temperature 
and  salinity  data  and  R  is  the  expected  error  from  the  01. 

The  observation  operator  //  links  the  state  vector  space  with  the  observation  space  and 
contains  the  individual  sound  speed  forecasts  of  the  NCOM  models. 

Similarly  to  the  LINREG  method,  the  SE  Kalman  filter  based  method  (KALMAN)  ean 
be  applied  to  models  resulting  from  an  EOF  regularization  and  is  denoted 
KALMAN  NORMEOF  hereafter.  This  approach  again  is  expected  to  increase 
generalization  skills  of  the  SE. 

In  oceanography,  usually,  the  state  vector  contains  hundred  of  thousands  of  grid  points, 
so  that  low-rank  approximations  of  the  Kalman  iteration  must  be  adopted  (e.g.  Pham  et 
al.,  1998).  Here,  the  state  veetor  is  very  small  (i.e.  the  three  weights,  two  on  the 
numerical  models  and  one  on  the  independent  term,  which  may  be  reduced  down  to  2  or 
1  after  the  PCA  regularization),  and  hence  the  full  Kalman  is  applied  at  every  grid  point. 
The  only  assumption  required  for  this  method  is  that  of  a  linear  model  transforming  the 
weights  in  time  and  Normal  weight  distributions. 

The  methods  described  above  have  been  tested  on  0-24hrs  forecast  for  3  consecutive  days 
with  a  running  window  from  29  April-  IMay  with  3  hrs  steps.  The  training  period  starts 
on  April  16  and  ends  on  29  April  OOhOO.  Figure  5  shows  resulting  sound  speed  deviations 
from  the  01  analyzed  fields  for  days  29  April  2007  12h00.  One  notes  that  both 
NCOM_COARSE  and  NCOM_FINE  are  rather  well  tuned  in  the  middle  of  the  water 
column,  except  that  the  fine  resolution  NCOM  shows  some  larger  ‘spotty’  errors. 
However  strong  biases  ean  be  identified  in  the  upper  15m  at  the  surface  and  the  bottom 
10m,  mainly  because  only  few  CTD  reached  the  bottom  layers.  The  ENSMEAN  is 
strongly  biased  as  well.  The  UNBIASED  ENSMEAN  seems  to  eorreet  most  of  the 
surface  and  bottom  bias  but  has  larger  errors  in  the  15-60m  range.  The  LINREG  and 


KALMAN  forecasts  correct  for  most  of  the  errors,  the  larger  discrepancies  at  the  surface 
having  been  decreased  significantly.  Qualitatively,  the  EOF  regularization  benefit  for 
KALMAN  is  not  directly  apparent  from  the  figure.  Results  for  1  May  12h00  arc 
qualitatively  similar  except  that  here,  the  benefit  of  the  regularization  on  the  KALMAN 
method  is  clearly  apparent. 
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Fig  5.  24hrs  sound  speed  (m/s)  forecast  errors  on  section  A-B  on  29  April  12h00  for 
the  various  models.  From  left  to  right  and  top  to  bottom:  NCOM_COARSE, 
NCOM_FINE,  ENSMEAN,  UNBIASED  ENSMEAN,  LINREG_NORM, 
LINREG_NORMEOF,  KALMAN_NORM  and  KALMAN_NORMEOF 
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Fig.  6  Same  as  figure  5  but  for  1  May  12h00. 

The  relative  weighting  on  the  2  numerical  models  and  the  independent  terms,  and  on  the 
EOF  compressed  models,  were  examined  carefully  and  show  a  strong  variability,  both  in 
spaee  (vertical  and  horizontal)  and  time,  with  values  usually  in  the  range  [-2  2]. 

The  weights  on  NCOM_COARSE  and  NCOM  FINE  for  the  LINREG  method  (Fig.  7) 
suggest  that  the  coarse  NCOM  model  is  usually  more  reliable  and  contribute  significantly 
to  the  SE  skill:  weights  are  usually  closer  to  one  over  the  whole  section  AB.  Contribution 
from  NCOM  FINE  are  smaller  and  sometimes  negative,  indicating  than  this  model  may 
be  out-of-phase  at  specific  locations.  Note  also  some  higher  magnitudes  of  the  weights, 
up  to  2,  below  100m  around  latitude  42.6 1°N,  illustrating  the  poor  correlations  between 
numerical  models  and  data  at  this  particular  location. 
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Fig  7  Weights  on  NCOM_COARSE  (top)  and  NCOMFINE  (bottom)  for  the 
LINREG  method  on  section  AB  versus  depth  (m). 

Interpretation  of  these  weights  in  physical  terms,  especially  for  the  KALMAN-based 
methods,  remains  difficult  as  illustrated  in  Fig  8.  The  weights  on  NCOM  COARSE 
remain  usually  stronger  and  contribute  more  to  the  SE  skills,  especially  on  the  upper  and 
lower  10m  of  the  water  column  in  the  first  half  of  the  learning  period.  One  also  notes  the 
important  role  of  the  independent  term  which  evolves  significantly  in  time.  This  may  be 
an  indication  of  the  strong  variability  of  skills  of  the  underlying  numerical  models. 


NCOM  COARSE 


19  Apr  21  Apr  23  Apr  25  Apr  27  Apr  29  Apr  19  Apr  21  Apr  23  Apr  25  Apr  27  Apr  29  Apr 


Fig.  8  Evolution  of  weights  on  NCOM  COARSE  (top),  NCOM  FINE  (middle)  and  the 
independent  term  (bottom)  over  the  learning  period  versus  depth  (m)  for  the 
KALMAN  NORM  method  for  a  profile  at  42.64°N  (middle  of  seetion  AB). 

Error  statistics  between  the  various  models  and  the  01  analysis  have  been  computed  for 
the  whole  29  April- 1  May  forecast  period,  and  arc  shown  in  Fig.  9.  As  expeeted,  root 
mean  square  (RMS)  errors  for  NCOM  COARSE,  NCOM  FINE  and  ENSMEAN  arc 
quite  large,  mainly  due  the  model  biases  at  the  top  15m  and  the  bottom  10m  of  the  water 
eolumns.  NCOMFINE  errors  are  also  larger  than  NCOM  COARSE.  The 
UNBIASED_ENSMEAN  method  corrects  for  a  large  portion  of  the  bias  and  illustrates 
the  benefit  of  the  method.  Further  methods  decrease  the  RMS  error  even  further,  the 
optimal  approach  being  the  KALMAN_NORMEOF.  Hence,  by  adapting  the  weights 
dynamically  on  EOF  compressed  models,  the  generalization  skills  of  the  SE  have  been 
maximized.  NCOMCOARSE  shows  very  good  correlation  to  in-situ  analyzed  data 
contrary  to  NCOM  FINE.  This  illustrates  the  fact  that  most  of  the  error  in 
NCOM  COARSE  is  due  to  the  bias,  whilst  NCOM  FINE,  although  better  suited  for 
eddies  and  fine-seale  processes,  is  not  able  to  eorrcct  for  the  phase  of  such  features.  It 
should  be  stressed  at  this  point  that  the  NCOM  runs  were  not  assimilating  any  local 
profile  data  whieh  are  used  here  as  independent  validation  set  by  comparison  with 
numerical  and  SE  numerical  foreeasts.  Note  however  that  the  analysis  is  smoothing  the 
observations  to  scales  that  arc  more  comparable  to  NCOM  COARSE  (2km-3hrs),  hence 
not  resolving  the  fine  structures  eventually  simulated  by  the  NCOM  FINE  runs. 


Validation  against  in-situ  sparse  data  was  tested  but  conclusions  were  less  conclusive  due 
to  the  scarce  and  inhomogeneous  data  distribution. 

The  bias  of  all  SE  methods  using  observations  is  lower  than  the  numerical  models  and 
their  ensemble  average.  The  standard  deviation  of  the  forecasts  also  suggests  that  SE 
methods,  and  KALMAN  NORMEOF  in  particular,  result  in  better  representation  of 
energy  spectrum  of  the  predictions.  The  synthetic  Taylor  skill  score  (Taylor,  2001)  is  a 
function  of  both  signal  energy  and  correlation  and  shows  the  improvements  in  skill 
obtained  by  using  progressively  more  complex  methods,  starting  from  the  numerical 
models  to  the  dynamic  weight  adjustment  with  regularization. 
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Fig.  9.  Error  statistics  for  24  hrs  sound  speed  (m/s)  forecast  from  29  April  to  1  May 
2007  for  the  8  prediction  methods.  From  left  to  right  and  top  to  bottom:  RMS  error, 
correlation,  bias,  energy  difference  and  Taylor  skill.  Bars  from  left  to  right 
correspond  respectively  to  NCOM_COARSE,  NCOM_FINE,  ENSMEAN, 
UNBIASED_ENSMEAN,  LINREG_NORM,  LINREG_NORMEOF, 

K ALM AN_N ORM  and  KALMAN_NORMEOF. 

3.  Ocean-Acoustics  sensitivity  analysis 

Ocean  modeling  outputs  can  be  used  for  a  number  of  applications,  among  which  the 
assessment  of  acoustic  propagation  constitutes  an  important  matter  of  interest.  For 
shallow  water  environments  in  particular,  it  is  well  known  that  the  environment  strongly 
influences  the  propagation  of  sound.  A  proper  description  of  the  environment  is  thus 


required  comprising  the  water-eolumn  aeoustie  properties  (sound  speed  profile  versus 
time  and  range,  sea-surfaec  roughness)  and  the  seabottom  geoacoustic  properties.  It  is 
well  established  that  for  shallow  waters,  a  precise  assessment  of  the  range-dependent 
seafloor  and  subseafloor  properties  is  essential,  whilst  the  acoustical  field  is  less  sensitive 
to  the  water  column  properties.  As  a  result,  the  oeean-aeousties  community  has  spent 
mueh  efforts  in  the  last  decade  to  properly  assess  the  seabottom  properties.  Recently, 
emphasis  was  put  on  comparing  the  relative  effects  of  uncertain  seabottom  description 
versus  uncertain  water  column  description.  Though,  no  general  statement  about  the 
importance  of  one  environmental  parameter  over  the  other  ean  be  made  (the  basic 
eonelusion  being  that  it  is  generally  ease-dependent),  it  was  nevertheless  concluded  that 
effective  modeling  of  the  water  eolumn  is  an  important  component  of  acoustic  REA  even 
though  the  seabottom  geoacoustic  properties  prevails  in  most  shallow  water 
environments.  Interestingly,  such  statements  had  already  been  made  previously  for 
studies  devoted  to  the  assessment  of  the  seabottom  properties  in  which  gcoacoustic 
inversion  methods  were  reported  to  have  performance  limitations  due  to  inadequate 
description  of  the  water  column  and  in  particular,  the  range  dependence  of  the  sound 
speed  profile  and  its  time  variability  during  the  measurements.  In  the  light  of  these 
observations,  the  ocean  modeling  approaches  presented  in  the  previous  sections 
constitute  an  interesting  test  ease  in  order  to  investigate  the  impaet  of  ocean  modeling 
upon  acoustic  propagation  assessment.  More  precisely,  the  AB  transect  of  the 
MREA/BP07  sea  trial  is  characterized  by  a  range-  and  time-dependent  sound  speed 
profile,  diversely  modeled  as  illustrated  in  the  previous  sections.  Transposing  the  ocean 
predictions  in  acoustic  predictions,  i.e.  mapping  ocean  modeling  into  acoustic  modeling, 
provides  a  way  to  estimate  how  mueh  aeousties  foreeasts  ean  be  affected  by  oeean 
predictions.  Subsequently,  the  idea  is  to  adopt  the  end-user  viewpoint  of  an  acoustician  to 
estimate  the  neeessary  and  sufficient  degree  of  sophistication  of  oeean  modeling  required 
for  aeoustie  modeling.  In  other  words,  is  there  an  interest  to  model  precisely  the  ocean 
conditions  when  an  approximate  description  of  the  water  column  is  already  available? 
Moreover,  another  viewpoint  that  is  worth  being  investigated  with  the  framework  of  the 
BP07  sea  trial  is  closely  linked  to  one  major  objective  of  the  experiment,  i.e.  the 
assessment  of  the  seabottom  properties  by  means  of  fused  scismie  imaging  and 
gcoacoustic  inversion  techniques.  Considering  that  the  sampling  of  the  seawater 
properties  is  relatively  searee  both  in  time  and  in  geographical  extension,  one  may 
wonder  whether  available  4D  ocean  modeling  fields  could  be  directly  used  as  inputs  to 
the  gcoaeoustic  inversion  algorithm. 

To  analyze  the  two  viewpoints,  a  sensitivity  study  based  on  intensive  aeoustie 
propagation  simulation  was  conducted.  For  the  following  simulations,  the  sea  surface  was 
considered  to  be  perfectly  flat  and  the  seabottom  properties  were  kept  constant  for  every 
acoustic  run.  A  range-dependent  bathymetry,  stratigraphic  and  aeoustie  model  of  the 
sediment  along  the  AB  transeet  was  taken  from  previous  seabottom  characterization 
studies  along  the  same  transect  (Hcrmand  and  Gerstoft,  1996;  Hermand,  1999)  and  an 
additional  geophysical  survey  of  the  transect  with  a  multibeam-echo  sounder  and  a 
seismic  subbottom  profiling  system  (Fig.  10).  As  in  Hermand  and  Gerstoft  (1996)  and 
Hermand  (1999),  the  stratigraphy  of  the  seabottom  was  simplified  in  order  to  keep  the 
most  important  features  that  have  an  impaet  at  aeoustie  frequencies  from  few  hundreds 


Hertz  (300)  to  few  kHz  (1.8).  The  seabottom  was  modeled  by  an  overlying  very  soft 
sediment  layer  of  clayed  unconsolidated  sediments  with  a  variable  thickness  (from  5  m  to 
9  m),  over  a  silty  harder  sediment.  The  geoaeoustie  parameters  of  those  layers  are 
illustrated  in  Fig.  11.  To  compute  the  complex  acoustic  pressure  fields,  the  range- 
dependent  propagation  model  based  on  the  parabolic  approximation  RAM  (Collins,  1989) 
was  employed.  The  complex  acoustic  pressure  fields  were  computed  for  a  series  of  9 
frequencies  from  300  Hz  to  1800  Hz  along  the  AB  transect,  and  for  each  available  time 
step  of  the  ocean  modeling  (24  steps  with  3-h  interval  from  April  29  to  May  1).  A  typical 
environmental  model  that  was  used  is  presented  in  Fig.  1 1. 
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Fig.  10.  Seismic  survey  of  the  AB  transect  using  a  seismic  towed  sub-bottom  profiler 


Fig.  11  Typical  scenario  for  the  acoustic  prediction  runs  along  the  AB  transect.  A 
range-dependent  sound  speed  profile  from  the  ocean-modeling  work  is  assumed. 
The  seabottom  is  composed  of  a  range-dependent  overlaying  clayed  layer  over  a 
silty  basement. 

Though  sensitivity  studies  are  quite  often  applied  in  oeean-aeousties  to  examine  the  effect 
of  various  environmental  parameters  on  the  acoustic  propagation  in  a  particular 
environment,  very  little  has  been  done  to  provide  a  quantitative  approach  to  that  issue. 
Typically,  the  sensitivity  analysis  is  often  earried  out  qualitatively  by  observing  the 
change  in  aeoustie  fields  that  results  from  an  environmental  change  (e.g.  Fig.  12).  In 


Dosso  et  al  (2007  and  references  herein),  quantitative  approaches  have  been  proposed 
based  on  normalized  bias  estimates  between  the  reference  transmission  losses  field  and 
the  perturbed  ones.  In  the  present  paper,  a  different  approach  based  on  the  cross- 
correlation  of  the  fields  is  applied.  In  Dosso  et  al  (2007),  it  is  shown  that  the  sensitivities 
are  range,  depth  and  frequency  dependent.  Here,  the  depth  dependency  is  not 
investigated.  Instead,  the  measure  of  sensitivity  consists  in  a  normalized  cross-correlation 
of  the  acoustic  pressure  fields  at  every  range  along  a  synthetic  vertical  line  array  (VLA) 
with  an  inter-element  spacing  of  1  m  for  a  set  of  frequencies.  The  VLA  coverage  is 
limited  to  the  water  column.  This  metric  is  related  to  the  classical  “Bartlett  processor” 
applied  within  the  gcoacoustic  inversion  framework  based  on  matched  field  processing 
(MFP)  (Tolstoy,  2003).  Two  normalized  cross-correlations  have  been  used  in  the 
analysis: 
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where  p  is  the  reference  acoustic  pressure  field  using  the  ocean  optimally  interpolated 
data  presented  in  the  previous  section,  q  is  the  pressure  field  with  the  ocean  forecast 
under  investigation,  Nz  the  number  of  hydrophones  of  the  synthetic  VLA,  p*  is  the 
conjugate  of  the  complex  acoustical  pressure  field  p.  As  this  point,  it  should  he  mentioned 
that  the  optimal  interpolation  field  is  considered  as  the  best  reference  water  column 
model.  It  is  remembered  though  that  it  is  limited  in  range  and  time  resolution. 

Both  measures  are  bounded  within  the  interval  [0  1],  1  corresponding  to  a  perfect  match 
of  the  compared  pressure  fields.  <(>i  investigates  the  correlation  of  the  vertical  spatial 
structure  of  the  acoustic  pressure  fields  in  phase  and  amplitude,  while,  <j>2  investigates 
only  the  vertical  spatial  structure  of  the  amplitude  of  the  acoustic  pressure  fields  (or  the 
transmission  loss). 

The  sensitivity  measures  defined  so  far  consider  changes  in  the  acoustic  field  at  a  fixed 
range  in  space.  However,  the  acoustic  field  perturbation  due  to  an  environmental 
perturbation  generally  includes  a  component  representing  a  spatial  shift  of  the  field  in 
addition  to  a  change  to  the  shifted  field  ( e.g .  Dosso  et  al,  2007).  By  not  taking  into 
account  that  field  shifting  effect,  our  approach  can  lead  to  over-estimated  impacts  of  the 
environment  upon  the  acoustics,  whilst  an  analysis  of  the  acoustic  fields  shows 
qualitatively  that  the  impact  may  be  lower.  To  avoid  that  effect,  a  numerical  spatial  field 
shifting  was  applied.  A  range-depth  window  of  preselected  size  is  defined  for  the 
reference  field,  centered  on  the  current  range.  The  tested  field  is  correlated  with  each  of 
the  combination  of  the  reference  field  inside  the  window.  Finally,  only  the  maximum  of 
the  obtained  results  is  retained.  For  the  study,  the  acoustic  field  was  computed  with  an 
elementary  cell-size  of  20m  x  lm.  The  reference  window  was  defined  +/-  100m  about  the 
current  range  and  the  vertical  shift  was  +/-  5m  about  the  current  depth  (i.e.  for  each  cell 


in  the  frequency  x  range  domain,  121  combinations  of  the  pressure  fields  were 
considered). 
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Fig.  12  Example  of  the  acoustical  pressure  fields  obtained  with  the  various  ocean 
modeling  outputs  and  the  optimally  interpolated  data  field. 


The  results  of  the  sensitivity  analysis  ean  be  considered  from  several  perspectives 
presented  in  the  following  subsections. 

3.1  Sensitivity  analysis  from  an  acoustic  propagation  perspective 


In  Fig.  13a, b,  the  resulting  <j>i  and  <f>2  fields  in  the  frequency  times  range  domain  are 
plotted  for  the  oeean-modeling  outputs  produced  for  1  May,  12:00.  Some  trends  ean  be 
elearly  identified.  They  obviously  overcome  the  sole  framework  of  super-ensemble  ocean 
modeling  and  ean  be  generalized  to  the  need  of  a  proper  description  of  range-dependent 
properties  of  the  water  column  as  illustrated  below: 

•  The  correlation  levels  are  frequency-  and  range-dependent. 

•  The  lower  the  frequency,  the  lower  the  sensitivity.  At  very  low  frequency 
(typieally  from  300  to  500  Hz),  little  effect  on  the  aeoustie  field  ean  be  observed, 
meaning  that  a  relatively  eoarse  description  of  the  water  eolumn  is  aeeeptable.  At 
higher  frequencies  this  is  no  longer  the  ease.  In  particular,  it  ean  be  seen  that 
aeeurate  aeoustieal  predictions  are  limited  to  the  very  elose  neighborhood  of  the 
sound  souree  at  those  frequencies.  This  is  a  well-known  effeet  related  to  the  faet 
that  the  higher  the  frequency,  the  lower  the  aeoustie  wavelength,  the  higher  the 
impaet  of  the  small  scale  environmental  features. 

•  The  greater  the  distance,  the  higher  the  sensitivity.  This  result  is  quite  trivial 
sinee,  the  aeoustie  propagation  at  long  ranges  is  directly  linked  to  the  past 
propagation.  Differences  tend  to  accumulate  with  increasing  range. 


Interestingly,  it  should  be  noticed  the  stronger  impact  of  the  environment  when  the  full 
complex  pressure  field  (amplitude  and  phase)  is  kept,  rather  than  when  only  the 
amplitude  is  considered.  Though  the  same  range-frequency  dependency  can  be  observed 
on  both  types  of  fields,  the  correlation  levels  arc  much  lower  when  the  phase  information 
is  kept,  which  demonstrates  how  much  this  information  is  sensitive  to  the  description  of 
the  environment.  As  soon  as  the  correlation  level  is  frequency  dependent,  it  is 
emphasized  that  an  accurate  range  dependent  description  of  the  water  column  properties 
needs  to  be  fed  to  the  acoustic  models  involved  in  applications  based  on  full  field 
estimates,  e.g.  matched-field  processing  (MFP),  and  model-based  matched  filter  (MBMF) 
processing,  the  latter  involving  the  prediction  of  the  band-limited  impulse  response  of  the 
medium.  The  importance  of  such  a  need  increases  with  frequency  and  range. 

In  Figs.  14a,b,  a  multi-frequency  correlation  measure  was  deduced  from  the  2D 
frequency  times  range  intercorrclation  matrices.  This  is  straightforwardly  obtained  by 
computing,  for  each  range  and  each  model,  the  mean  of  the  <f>i  and  fc  outputs  over 
frequency,  which  leads  to  the  generalized  multi- frequency  Bartlett  processor  outputs  that 
are  routinely  in  acoustic  inversion  works  based  on  matched  field  processing  (MFP). 
Whatever  the  frequency  and  the  range,  multi-frequency  correlations  of  the  amplitude 
fields  remain  relatively  high  (over  0.85)  reinforcing  the  conclusion  that  for  these 
frequencies,  ocean  modeling  efforts  might  have  a  limited  scope  (no  measure  of  the  degree 
of  adequacy  of  the  water  column  model  is  discussed  here  though).  For  full  complex 
pressure  fields,  the  decay  of  the  inter-correlation  fields  is  much  higher  with  range.  Only 
the  two  KALMAN  filter-based  models  arc  shown  to  maintain  good  performances 
(correlation  values  higher  than  0.7  for  ranges  up  to  5  km),  whilst  the  remaining  models 
are  much  more  limited  in  range.  Here  also,  the  conclusions  arc  identical:  the  computation 
of  the  full  complex  pressure  field  requires  a  more  accurate  description  of  the  water 
column  conditions.. 


Fig.  13.  2D  sensitivities  of  (a)  the  complex-valued  acoustic  pressure  field  for  the 
eight  ocean  modeling  outputs  and  (b)  the  amplitude  of  the  acoustic  field  on  1  May 
2007,  12h00. 


Fig.  14  Multi-frequency  sensitivities  for  1  May  2007  at  12h00:  (a)  complex-valued 
acoustic  pressure  field  and  (b)  corresponding  amplitude. 


3.2  Sensitivity  analysis  from  an  ocean-modeling  perspective 

From  an  ocean  modeling  perspective,  the  sensitivity  analysis  allows  to  extract  interesting 
features  about  the  performances  of  the  different  environmental  forecasts. 


As  expected  from  the  evaluation  of  the  super-ensemble  ocean-modeling  skills  presented 
in  the  previous  sections,  better  matches  of  the  acoustic  fields  can  generally  be  obtained 
when  the  super-ensemble  outputs  are  employed  in  place  of  the  initial  NCOM  simulations. 
In  some  instances  though,  the  improvements  arc  quite  low  especially  for  the  less 
sophisticated  super-ensemble  techniques  (ENSMEAN,  UNBIASED  ENSMEAN).  The 
UNBIASED  ENSMEAN  technique  was  even  seen  to  lead  to  less  accurate  acoustic 
forecasts  for  the  very  first  predictions  on  29  April  (sec  Fig  15a). 


A  hierarchy  of  the  skills  of  super-ensemble  techniques  from  an  acoustic  modeling 
viewpoint  can  be  established  from  the  intercorrelation  fields.  The  improvement  given  by 
the  ENSMEAN  based  techniques  are  generally  lower  than  those  of  the  LINREG  ones, 
which  are  also  less  skilled  than  the  KALMAN-based  techniques.  More  interestingly,  it 
can  be  emphasized  that  the  KALMAN  filter  techniques  shows  rather  constant  consistency 
(both  in  range  and  frequency)  over  the  3  days  of  the  test  (sec  Fig.  14a  and  15a,b). 
Meanwhile,  the  two  other  types  of  methodologies  arc  far  less  robust  over  time.  For 
example,  the  UNBIASED  ENSMEAN  approach  starts  with  poorer  skills  than  any  other 
forecast  (sec  Fig.  15a),  and  about  the  middle  of  the  runs,  it  suddenly  improves  over  the 
other  methods,  except  the  KALMAN-based  forecasts  (sec  Fig.  14a).  The  reasons  for  such 
an  improvement  seem  to  be  linked  to  a  more  accurate  assessment  of  the  thcrmoeline  for 
this  method  after  several  runs. 


Though  the  most  sophisticated  methods  showed  better  matches  of  the  acoustic  fields  with 
respect  to  the  reference  water-column  field,  none  of  them  was  nevertheless  able  to  cope 
with  all  the  characteristics  of  this  reference  field,  with  a  higher  impact  at  long  ranges  and 
high  frequencies.  Moreover,  it  should  also  be  emphasized,  that  though  the  oceanographic 
survey  was  rather  important,  the  optimally  interpolated  data  still  present  some  limitations 
in  range  and  time  resolutions,  i.e.  2  km  and  1  day  respectively.  Therefore,  for  the  purpose 


of  acoustic  predictions,  there  is  a  need  for  improvement  of  the  ocean  forecasts.  Several 
ways  arc  possible  to  enhance  ocean  modeling  skills  for  the  specific  requirements  of 
acoustics.  A  denser  hydrographic  survey  and  a  better  assimilation  of  in  situ  ocean  data  at 
earlier  stages  of  the  ocean-modeling  phase  are  likely  to  provide  better  initial  forecasts. 
Anyway,  this  approach  can  be  constrained  by  the  number  of  measuring  platforms  and 
there  may  also  be  conflicts  between  the  measurement  strategics  for  ocean-modeling 
solely  and  those  dedicated  to  acoustics.  An  alternate  approach  which  could  be 
particularly  worthwhile  within  the  current  ocean-acoustics  framework  is  to  use  acoustic 
observations  to  resolve  range-dependent  oceanographic  fields  that  can  be  used  in  a 
second  step  to  improve,  within  a  feedback  loop  the  ocean  modeling  outputs.  The 
interesting  aspect  of  this  approach  is  that  acoustic  observations  lead  to  synoptic 
observations  of  the  ocean  that  presumably  should  allow  catching  the  oceanic  variability  at 
fine  scales  and  long  ranges.  Moreover,  these  measurements  are  more  closely  linked  to  the 
final  application.  The  most  relevant  features  of  the  environment  are  likely  to  be  caught 
through  the  acoustic  observation  and  thus  the  resulting  ocean  modeling  is  expected  to 
more  closely  serve  the  acoustic  needs.  A  companion  paper  investigates  the  feasibility  of 
such  a  range-dependent  tomography  scheme  (Carrierc  et  al,  this  issue)  and  discusses 
some  of  the  basis  for  such  a  general  ocean-acoustics  feedback  loop. 


( 
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Fig.  15  Multi-frequency  sensitivities  on  (a)  29  April  2007,  12h00  and  (b)  1  May  2007, 
OOhOO. 

3.3  Sensitivity  analysis  from  a  seabottom  geoacoustic  inversion  perspective 

The  use  of  4D  ocean  modeling  outputs  constitutes  an  interesting  alternative  for  the 
gcoacoustic  inversion  community.  Since  range  dependency  of  the  water  column  affects 
acoustic  propagation,  any  inversion  scheme  for  seabottom  acoustic  properties  has  to 
constrain  the  environmental  model  with  measured  or  predicted  sound  speed  profiles  or  to 
search  for  the  best  (ad-hoc)  sound  speed  field  in  searching  for  the  seabottom  properties 
(range-dependent  geoacoustic  inversion).  Another  viewpoint  from  the  inversion 
perspective  is  to  investigate  the  problem  with  a  “system-oriented”  approach,  i.e.  what  are 
the  impact  of  range  dependency  on  the  inversion  results  and  whether  we  can  cope  with 
water  column  range  dependency  when  performing  geoacoustic  inversion  with  a  range 


independent  model.  Since  range  dependency  of  the  water  eolumn  ean  limit  the 
performances  of  the  range-independent  approach  efforts  have  been  spent  on  developing 
the  range-dependent  approach.  However  the  resulting  methods  are  often  quite  time 
consuming  beeause  they  require  eomputationally-expensive  forward  aeoustic  propagation 
models.  Moreover,  they  often  require  good  a  priori  knowledge  of  the  range  dependency 
to  better  pose  the  inverse  problem.  To  simplify  the  inverse  problem,  it  can  thus  be  useful 
to  limit  the  inverse  problem  to  a  geometrical  set  up  that  limits  the  impact  of  the  range 
dependency.  Here  it  appears  that  at  close  ranges  i.c.,  from  a  few  hundreds  meters  to  a  few 
kilometers,  the  correlations  between  predieted  and  measured  fields  remain  high  over  the 
whole  water  column  (over  0.8-0.9),  whatever  the  water  column  conditions  arc.  This  gives 
confidence  that  experimental  configurations  involving  short  ranges  arc  most  likely  to  be 
weakly  impacted  by  range  dependency.  Such  an  approach  has  been  followed  during  the 
BP07  sea  trial  during  which  dedicated  seabottom  geoacoustic  inversion  runs  were 
designed  to  assess  the  seabottom  properties  in  parallel  to  the  ocean  modeling  efforts.  For 
these  runs,  a  sound  source  emitting  broadband-coded  signals  within  the  frequency  range 
300-1800  Hz  was  employed  and  the  signals  were  received  along  sparse  VLAs  at 
distances  not  exceeding  1-2  km.  Though  the  full  vertical  intcrcorrelation  measure  is  not 
completely  representative  of  the  VLA  design  employed  (sparse  arrays  limited  to  a  part  of 
the  water  column),  it  provides  good  confidence  that  such  a  setup  is  suited  for  gcoacoustic 
inversion  purposes  with  little  impaet  of  the  range  dependency  of  the  water  column  as  was 
recently  confirmed  by  the  first  BP07  geoacoustic  inversion  results  along  the  AB  transect 
(Hcrmand  and  Lc  Gac,  2008).  Hopefully,  this  is  expected  to  be  extended  to  the  small- 
seale  time  variability.  Further  work  shall  be  dedicated  to  investigate  that  assumption. 


4.  Conclusions  and  perspectives 

This  work  presents  an  innovative  approach  to  the  eouplcd  ocean-acoustic  problem. 
Despite  scarce  hydrographic  monitoring  which  is  obviously  a  strong  limitation  for 
aeeurate  and  reliable  oeean  foreeasting,  sophisticated  -  but  ehcap  -  methods  ean  be  used 
to  exploit  both  observational  and  prediction  information  streams.  We  have  demonstrated 
the  potential  benefit  of  combining  dynamically  various  observational  data  sourees  (i.e. 
CTD,  thermistor  strings  and  others)  and  traditional  oeean  prediction  into  a  super- 
ensemble  forecast  to  better  characterize  and  predict  ocean  and  acoustic  properties. 

In  a  companion  paper  (Carricrc  et  al,  this  issue),  tomography-derived  synoptic  data 
provides  a  promising  complement  (or  even  an  alternative  approach)  to  standard 
hydrographic  measurements.  In  particular,  range  dependent  temperature  or  sound  speed 
profile  obtained  by  acoustic  inversion  may  be  in  turn  assimilated  in  ocean  models  or 
multi-model  super-ensembles  (SE).  Subsequent  nowcast  and  forecast  can  then  further 
help  refining  the  inversion  proeess. 

Ongoing  efforts  on  SE  include  the  exploitation  of  remote  sensing  information  to  better 
understand  SE  error  eovariance  and  to  overcome  the  lack  of  in-situ  observations. 


The  prediction  of  the  SE  error  statistics  also  deserves  some  increased  efforts.  The 
uncertainties  may  be  exploited  in  the  inversion  process  because  they  help  better 
constraining  the  search  space.  These  uncertainties  may  also  be  used  in  adaptive  sampling 
efforts  to  optimize  the  observational  network  and  asset  allocation  by  focusing  on  areas 
where  errors  arc  expected  to  be  higher. 

This  multidisciplinary  ocean-acoustic  feedback  approach  offers  some  promising 
perspectives  in  the  context  of  Rapid  Environmental  Assessment. 
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